AlpineStreamInsectMicrobiome
JReceveur
Data overview
Sequencing
physeq## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2420 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 2420 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 2420 tips and 2400 internal nodes ]Metadata
Trtdata <- ddply(metadata, c("Sampling_station", "FFG","Taxon_name"), summarise,
N = length(id),
meanWeight = mean(Mass),
sd = sd(Mass),
se = sd / sqrt(N)
)
Trtdata## Sampling_station FFG Taxon_name N meanWeight sd
## 1 Ostana Filterer Hydropsyche 4 41.4500 33.707615
## 2 Ostana Predator Perla 2 190.5000 174.513954
## 3 Ostana Scraper Epeorus 5 22.7600 3.857849
## 4 Ostana Shredder Tipula_(Acutitipula) 3 534.9333 425.432889
## 5 PDR Predator Dictyogenus 4 83.9500 16.270730
## 6 PDR Scraper Epeorus 4 11.6000 4.009988
## 7 PDR Shredder Tipula_(Acutitipula) 4 87.7000 42.817987
## se
## 1 16.853808
## 2 123.400000
## 3 1.725283
## 4 245.623793
## 5 8.135365
## 6 2.004994
## 7 21.408993ggplot(metadata,aes(x=Sampling_station,y=Mass,fill=Sampling_station))+geom_boxplot()+ #Write in labels from posthoc
scale_fill_manual(values=cbPalette)+xlab("")+ylab("Mass (mg)")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG,scales = "free_y")#+theme(axis.text.x = element_text(angle = 45, hjust = 1))+#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
head(metadata)## id FieldID Taxon_name FFG Sampling_station Mass
## 1 II5 T1OS Tipula_(Acutitipula) Shredder Ostana 989.1
## 2 II6 T2OS Tipula_(Acutitipula) Shredder Ostana 470.0
## 3 II8 T4OS Tipula_(Acutitipula) Shredder Ostana 145.7
## 4 II13 P1OS Perla Predator Ostana 67.1
## 5 II14 P2OS Perla Predator Ostana 313.9
## 6 II15 H2OS Hydropsyche Filterer Ostana 7.2
## shannon faith_pd
## 1 7.086586 20.95169
## 2 7.430371 22.58281
## 3 7.383975 21.47860
## 4 6.952585 17.94331
## 5 6.200620 26.41718
## 6 6.153952 17.75811metadataStatsSubset<- subset(metadata,FFG!="Predator")
metadataStatsSubset<- subset(metadataStatsSubset,FFG!="Filterer")
metadataStatsSubset## id FieldID Taxon_name FFG Sampling_station Mass
## 1 II5 T1OS Tipula_(Acutitipula) Shredder Ostana 989.1
## 2 II6 T2OS Tipula_(Acutitipula) Shredder Ostana 470.0
## 3 II8 T4OS Tipula_(Acutitipula) Shredder Ostana 145.7
## 10 E1OS E1OS Epeorus Scraper Ostana 20.2
## 11 E2OS E2OS Epeorus Scraper Ostana 27.5
## 12 E3OS E3OS Epeorus Scraper Ostana 19.0
## 13 E4OS E4OS Epeorus Scraper Ostana 20.8
## 14 E5OS E5OS Epeorus Scraper Ostana 26.3
## 15 II1 T1P Tipula_(Acutitipula) Shredder PDR 131.8
## 16 II2 T2P Tipula_(Acutitipula) Shredder PDR 53.9
## 17 II3 T3P Tipula_(Acutitipula) Shredder PDR 48.2
## 18 II4 T4P Tipula_(Acutitipula) Shredder PDR 116.9
## 23 E1P E1P Epeorus Scraper PDR 11.0
## 24 E3P E3P Epeorus Scraper PDR 17.4
## 25 E4P E4P Epeorus Scraper PDR 9.6
## 26 E5P E5P Epeorus Scraper PDR 8.4
## shannon faith_pd
## 1 7.0865864 20.951686
## 2 7.4303709 22.582812
## 3 7.3839746 21.478599
## 10 4.9093727 11.372657
## 11 1.9619653 9.905822
## 12 3.0828744 9.146140
## 13 4.8326128 7.641727
## 14 3.0611306 7.197870
## 15 7.2937707 21.728816
## 16 6.9992923 20.798101
## 17 6.0582633 11.170252
## 18 7.1744810 21.997409
## 23 1.6024907 4.266886
## 24 4.5283738 5.468116
## 25 0.7750815 4.605169
## 26 0.7526593 4.586619print("Mass (mg)~ Sampling_station*FFG, method=ANOVA")## [1] "Mass (mg)~ Sampling_station*FFG, method=ANOVA"av=aov( Mass~ Sampling_station*FFG, data=metadataStatsSubset)
summary(av)## Df Sum Sq Mean Sq F value Pr(>F)
## Sampling_station 1 109131 109131 3.563 0.08351 .
## FFG 1 319410 319410 10.427 0.00723 **
## Sampling_station:FFG 1 184026 184026 6.007 0.03054 *
## Residuals 12 367594 30633
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("(Since the interaction between sampling station and FFG was signinficant below are t-tests split out by FFG/n Mass (mg) ~Sampling-station (grouped by FFG, FDR correction)")## [1] "(Since the interaction between sampling station and FFG was signinficant below are t-tests split out by FFG/n Mass (mg) ~Sampling-station (grouped by FFG, FDR correction)"Scrapers<-subset(metadataStatsSubset, FFG=="Scraper")
t.test( Mass~ Sampling_station, data=Scrapers)##
## Welch Two Sample t-test
##
## data: Mass by Sampling_station
## t = 4.2191, df = 6.4396, p-value = 0.004753
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.79325 17.52675
## sample estimates:
## mean in group Ostana mean in group PDR
## 22.76 11.60Shredders<-subset(metadataStatsSubset, FFG=="Shredder")
t.test( Mass~ Sampling_station, data=Shredders)##
## Welch Two Sample t-test
##
## data: Mass by Sampling_station
## t = 1.8139, df = 2.0304, p-value = 0.2095
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -598.5184 1492.9851
## sample estimates:
## mean in group Ostana mean in group PDR
## 534.9333 87.7000#compare_means(Mass_.mg. ~ Sampling_station,group.by = "FFG", data = metadata, p.adjust.method = "fdr",method="t.test")Alpha Diversity
Observed Species
plot_richness(physeq, x="FFG",color="FFG", measures=c("Observed"))+geom_boxplot(aes(x=FFG, y=value, color=FFG), alpha=0.05)+ylab("Observed Species")+geom_point(size=3)+guides(fill=FALSE)
Shannon Diversity
av=kruskal.test( shannon~ FFG, data=metadataStatsSubset) #fig.width = 8, fig.height = 8 for line above to change size
av2=kruskal.test( shannon~ Sampling_station, data=metadataStatsSubset)
av##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 11.118, df = 1, p-value = 0.0008551av2##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 0.54044, df = 1, p-value = 0.4622posthoc<-compare_means(shannon ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(metadataStatsSubset,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
stats2<-summarySE(metadata,measurevar="shannon",groupvars=c("FFG","Sampling_station",na.rm=TRUE))
posthoc2<-compare_means(shannon ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr",group.by = "Sampling_station")
ggplot(metadata,aes(x=FFG,y=shannon,fill=FFG))+geom_boxplot() +#Write in labels from posthoc+geom_col()+geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black")+
scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~Sampling_station)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
ggplot(metadata,aes(x=Sampling_station,y=shannon,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
av##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 11.118, df = 1, p-value = 0.0008551posthoc## # A tibble: 1 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 shannon Scraper Shredder 0.00103 0.001 0.001 ** WilcoxonLetters## Scraper Shredder
## "a" "b"stats## FFG na.rm N shannon sd se ci
## 1 Scraper NA 9 2.834062 1.6633640 0.5544547 1.2785748
## 2 Shredder TRUE 7 7.060963 0.4686418 0.1771300 0.4334214posthoc2## # A tibble: 2 x 9
## Sampling_station .y. group1 group2 p p.adj p.format p.signif
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
## 1 Ostana shan~ Scrap~ Shred~ 0.0369 0.037 0.037 *
## 2 PDR shan~ Scrap~ Shred~ 0.0304 0.037 0.030 *
## # ... with 1 more variable: method <chr>stats2## FFG Sampling_station na.rm N shannon sd se ci
## 1 Filterer Ostana NA 4 6.502903 0.3062089 0.1531045 0.4872467
## 2 Predator Ostana NA 2 6.576603 0.5317195 0.3759825 4.7773101
## 3 Predator PDR NA 4 6.164048 1.3841309 0.6920654 2.2024611
## 4 Scraper Ostana NA 5 3.569591 1.2718220 0.5687761 1.5791756
## 5 Scraper PDR NA 4 1.914651 1.7867880 0.8933940 2.8431785
## 6 Shredder Ostana TRUE 3 7.300311 0.1865387 0.1076982 0.4633878
## 7 Shredder PDR NA 4 6.881452 0.5619605 0.2809802 0.8942045av2##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 0.54044, df = 1, p-value = 0.4622AlphaDivSplitByFFG
print("Scrapers")## [1] "Scrapers"Scrapers<-subset(metadataStatsSubset, FFG=="Scraper")
shapiro.test(Scrapers$shannon)##
## Shapiro-Wilk normality test
##
## data: Scrapers$shannon
## W = 0.8943, p-value = 0.2208shapiro.test(Scrapers$faith_pd)##
## Shapiro-Wilk normality test
##
## data: Scrapers$faith_pd
## W = 0.91387, p-value = 0.3439t.test( shannon~ Sampling_station, data=Scrapers)##
## Welch Two Sample t-test
##
## data: shannon by Sampling_station
## t = 1.5626, df = 5.2748, p-value = 0.1759
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.025406 4.335285
## sample estimates:
## mean in group Ostana mean in group PDR
## 3.569591 1.914651t.test( faith_pd~ Sampling_station, data=Scrapers)##
## Welch Two Sample t-test
##
## data: faith_pd by Sampling_station
## t = 5.3855, df = 4.8851, p-value = 0.003194
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.243935 6.398356
## sample estimates:
## mean in group Ostana mean in group PDR
## 9.052843 4.731698kruskal.test( shannon~ Sampling_station, data=Scrapers)##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 2.94, df = 1, p-value = 0.08641kruskal.test( faith_pd~ Sampling_station, data=Scrapers)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 6, df = 1, p-value = 0.01431print("Shredders")## [1] "Shredders"Shredders<-subset(metadataStatsSubset, FFG=="Shredder")
t.test(shannon~Sampling_station,data=Shredders)##
## Welch Two Sample t-test
##
## data: shannon by Sampling_station
## t = 1.392, df = 3.8225, p-value = 0.2394
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.4321367 1.2698544
## sample estimates:
## mean in group Ostana mean in group PDR
## 7.300311 6.881452t.test( faith_pd~ Sampling_station, data=Shredders)##
## Welch Two Sample t-test
##
## data: faith_pd by Sampling_station
## t = 1.0402, df = 3.2033, p-value = 0.3703
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.364469 10.859245
## sample estimates:
## mean in group Ostana mean in group PDR
## 21.67103 18.92364kruskal.test(shannon~Sampling_station,data=Shredders)##
## Kruskal-Wallis rank sum test
##
## data: shannon by Sampling_station
## Kruskal-Wallis chi-squared = 2, df = 1, p-value = 0.1573kruskal.test( faith_pd~ Sampling_station, data=Shredders)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.5, df = 1, p-value = 0.4795AlphaDivSplitByStation
print("Ostana")## [1] "Ostana"Ostana<-subset(metadata, Sampling_station=="Ostana")
OstanaStatAlpha<-subset(Ostana,FFG!="Predator")
shapiro.test(OstanaStatAlpha$shannon)##
## Shapiro-Wilk normality test
##
## data: OstanaStatAlpha$shannon
## W = 0.87409, p-value = 0.07366shapiro.test(OstanaStatAlpha$faith_pd)##
## Shapiro-Wilk normality test
##
## data: OstanaStatAlpha$faith_pd
## W = 0.8379, p-value = 0.02611kruskal.test( shannon~ FFG, data=OstanaStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 9.6923, df = 2, p-value = 0.007859compare_means(shannon ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 shannon Filterer Scraper 0.0200 0.052 0.020 * Wilcoxon
## 2 shannon Filterer Shredder 0.0518 0.052 0.052 ns Wilcoxon
## 3 shannon Scraper Shredder 0.0369 0.052 0.037 * Wilcoxonkruskal.test( faith_pd~ FFG, data=OstanaStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 8.2564, df = 2, p-value = 0.01611compare_means(faith_pd ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 faith_pd Filterer Scraper 0.0200 0.055 0.020 * Wilcoxon
## 2 faith_pd Filterer Shredder 0.596 0.6 0.596 ns Wilcoxon
## 3 faith_pd Scraper Shredder 0.0369 0.055 0.037 * Wilcoxonprint("PDR")## [1] "PDR"PDRStatAlpha<-subset(metadata, Sampling_station!="Ostana")
kruskal.test( shannon~ FFG, data=PDRStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: shannon by FFG
## Kruskal-Wallis chi-squared = 7.7308, df = 2, p-value = 0.02095compare_means(shannon ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 shannon Predator Scraper 0.0304 0.046 0.03 * Wilcoxon
## 2 shannon Predator Shredder 0.470 0.47 0.47 ns Wilcoxon
## 3 shannon Scraper Shredder 0.0304 0.046 0.03 * Wilcoxonkruskal.test( faith_pd~ FFG, data=PDRStatAlpha)##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 7.5385, df = 2, p-value = 0.02307compare_means(faith_pd ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")## # A tibble: 3 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 faith_pd Predator Scraper 0.0304 0.046 0.03 * Wilcoxon
## 2 faith_pd Predator Shredder 0.665 0.67 0.67 ns Wilcoxon
## 3 faith_pd Scraper Shredder 0.0304 0.046 0.03 * Wilcoxonposthoc<-compare_means(shannon ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(OstanaStatAlpha,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity (Ostana)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
posthoc<-compare_means(shannon ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(PDRStatAlpha,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity (PDR)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
posthoc<-compare_means(faith_pd ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(OstanaStatAlpha,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=faith_pd,fill=FFG))+geom_text(aes(x=FFG, y=faith_pd+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=faith_pd-se,ymax=faith_pd+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD (Ostana)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
posthoc<-compare_means(faith_pd ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(PDRStatAlpha,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=faith_pd,fill=FFG))+geom_text(aes(x=FFG, y=faith_pd+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+ geom_errorbar(aes(ymin=faith_pd-se,ymax=faith_pd+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD (PDR)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))
Faith’s PD
av=kruskal.test( faith_pd~ FFG, data=metadataStatsSubset)
av2=kruskal.test( faith_pd~ Sampling_station, data=metadataStatsSubset)
av##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 10.423, df = 1, p-value = 0.001245av2##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.89338, df = 1, p-value = 0.3446posthoc<-compare_means(faith_pd ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr")
Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
stats<-summarySE(metadata,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(metadata,aes(x=FFG,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD")+labs(fill="FFG")#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))geom_text(aes(x=FFG, y=faith_pd+se+2,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ 
stats2<-summarySE(metadata,measurevar="faith_pd",groupvars=c("FFG","Sampling_station",na.rm=TRUE))
posthoc2<-compare_means(faith_pd ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr",group.by = "Sampling_station")
ggplot(metadata,aes(x=FFG,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith PD")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~Sampling_station,scales = "free_x")+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
ggplot(metadata,aes(x=Sampling_station,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 
av##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by FFG
## Kruskal-Wallis chi-squared = 10.423, df = 1, p-value = 0.001245posthoc## # A tibble: 1 x 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 faith_pd Scraper Shredder 0.00150 0.0015 0.0015 ** WilcoxonLetters## Scraper Shredder
## "a" "b"stats## FFG na.rm N faith_pd sd se ci
## 1 Filterer NA 4 21.911376 2.894769 1.4473845 4.606223
## 2 Predator NA 6 21.557631 5.684157 2.3205474 5.965157
## 3 Scraper NA 9 7.132334 2.594210 0.8647367 1.994086
## 4 Shredder TRUE 7 20.101096 3.984788 1.5061083 3.685314posthoc2## # A tibble: 2 x 9
## Sampling_station .y. group1 group2 p p.adj p.format p.signif
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
## 1 Ostana fait~ Scrap~ Shred~ 0.0369 0.037 0.037 *
## 2 PDR fait~ Scrap~ Shred~ 0.0304 0.037 0.030 *
## # ... with 1 more variable: method <chr>av2##
## Kruskal-Wallis rank sum test
##
## data: faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.89338, df = 1, p-value = 0.3446Taxa Plots
GPr = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
GPr = filter_taxa(GPr, function(x) mean(x) > 1e-5, TRUE)
PhylumAll=tax_glom(GPr, "Phylum")
PhylumLevel = filter_taxa(PhylumAll, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower tha 0.1%
FamilyAll=tax_glom(GPr,"Family")
FamilyLevel = filter_taxa(FamilyAll, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%
GenusAll=tax_glom(GPr,"Genus")
GenusLevel = filter_taxa(GenusAll, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%Phylum Level
The top 7 phyla across all samples are given in the table below
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
df2 <- psmelt(PhylumAll)
df2$Abundance<-df2$Abundance*100
PhylumAll<-ddply(df2, c("Phylum"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
PhylumSorted<-PhylumAll[order(-PhylumAll$mean),]
PhylumSorted[1:7,]## Phylum N mean sd se
## 20 Proteobacteria 26 51.896154 21.817321 4.2787287
## 5 Bacteroidetes 26 17.466667 14.180099 2.7809462
## 13 Firmicutes 26 13.275641 17.088823 3.3513939
## 23 Tenericutes 26 10.424359 19.936032 3.9097776
## 26 Verrucomicrobia 26 1.524359 1.846861 0.3621992
## 19 Planctomycetes 26 1.494872 1.825686 0.3580464
## 3 Actinobacteria 26 1.121795 1.100947 0.2159134PhylumAll2<-ddply(df2, c("Phylum","FFG"), summarise, #For Rel abu tabs
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
PhylumAll2## Phylum FFG N mean sd se
## 1 [Thermi] Filterer 4 0.000000000 0.00000000 0.000000000
## 2 [Thermi] Predator 6 0.011111111 0.01721326 0.007027284
## 3 [Thermi] Scraper 9 0.000000000 0.00000000 0.000000000
## 4 [Thermi] Shredder 7 0.080952381 0.07663560 0.028965536
## 5 Acidobacteria Filterer 4 0.166666667 0.20905431 0.104527154
## 6 Acidobacteria Predator 6 0.472222222 0.37677973 0.153819680
## 7 Acidobacteria Scraper 9 0.111111111 0.14043583 0.046811943
## 8 Acidobacteria Shredder 7 0.376190476 0.28396289 0.107327883
## 9 Actinobacteria Filterer 4 0.383333333 0.48112522 0.240562612
## 10 Actinobacteria Predator 6 1.311111111 0.38624498 0.157683853
## 11 Actinobacteria Scraper 9 1.225925926 1.28914030 0.429713434
## 12 Actinobacteria Shredder 7 1.247619048 1.46980184 0.555532879
## 13 Armatimonadetes Filterer 4 0.033333333 0.04714045 0.023570226
## 14 Armatimonadetes Predator 6 0.255555556 0.13109228 0.053518198
## 15 Armatimonadetes Scraper 9 0.000000000 0.00000000 0.000000000
## 16 Armatimonadetes Shredder 7 0.104761905 0.11454996 0.043295815
## 17 Bacteroidetes Filterer 4 22.616666667 3.76696164 1.883480820
## 18 Bacteroidetes Predator 6 32.700000000 10.71460063 4.374217390
## 19 Bacteroidetes Scraper 9 0.844444444 0.65021364 0.216737880
## 20 Bacteroidetes Shredder 7 22.838095238 6.05320327 2.287895784
## 21 Chlamydiae Filterer 4 0.000000000 0.00000000 0.000000000
## 22 Chlamydiae Predator 6 0.000000000 0.00000000 0.000000000
## 23 Chlamydiae Scraper 9 0.018518519 0.03767961 0.012559870
## 24 Chlamydiae Shredder 7 0.000000000 0.00000000 0.000000000
## 25 Chlorobi Filterer 4 0.025000000 0.05000000 0.025000000
## 26 Chlorobi Predator 6 0.083333333 0.10697871 0.043673876
## 27 Chlorobi Scraper 9 0.000000000 0.00000000 0.000000000
## 28 Chlorobi Shredder 7 0.004761905 0.01259882 0.004761905
## 29 Chloroflexi Filterer 4 0.058333333 0.09574271 0.047871355
## 30 Chloroflexi Predator 6 0.066666667 0.06324555 0.025819889
## 31 Chloroflexi Scraper 9 0.003703704 0.01111111 0.003703704
## 32 Chloroflexi Shredder 7 0.000000000 0.00000000 0.000000000
## 33 Cyanobacteria Filterer 4 0.125000000 0.08333333 0.041666667
## 34 Cyanobacteria Predator 6 1.172222222 0.97920981 0.399760731
## 35 Cyanobacteria Scraper 9 0.162962963 0.29927897 0.099759656
## 36 Cyanobacteria Shredder 7 0.328571429 0.35194095 0.133021177
## 37 Deferribacteres Filterer 4 2.766666667 2.05065482 1.025327409
## 38 Deferribacteres Predator 6 0.016666667 0.04082483 0.016666667
## 39 Deferribacteres Scraper 9 0.000000000 0.00000000 0.000000000
## 40 Deferribacteres Shredder 7 0.814285714 0.57214902 0.216252002
## 41 Elusimicrobia Filterer 4 0.008333333 0.01666667 0.008333333
## 42 Elusimicrobia Predator 6 0.000000000 0.00000000 0.000000000
## 43 Elusimicrobia Scraper 9 0.000000000 0.00000000 0.000000000
## 44 Elusimicrobia Shredder 7 0.000000000 0.00000000 0.000000000
## 45 Fibrobacteres Filterer 4 0.000000000 0.00000000 0.000000000
## 46 Fibrobacteres Predator 6 0.005555556 0.01360828 0.005555556
## 47 Fibrobacteres Scraper 9 0.000000000 0.00000000 0.000000000
## 48 Fibrobacteres Shredder 7 0.000000000 0.00000000 0.000000000
## 49 Firmicutes Filterer 4 20.025000000 13.00792849 6.503964247
## 50 Firmicutes Predator 6 0.255555556 0.43851602 0.179023414
## 51 Firmicutes Scraper 9 0.788888889 1.23288280 0.410960934
## 52 Firmicutes Shredder 7 36.633333333 7.72724856 2.920625431
## 53 Fusobacteria Filterer 4 0.000000000 0.00000000 0.000000000
## 54 Fusobacteria Predator 6 0.000000000 0.00000000 0.000000000
## 55 Fusobacteria Scraper 9 0.003703704 0.01111111 0.003703704
## 56 Fusobacteria Shredder 7 0.000000000 0.00000000 0.000000000
## 57 Gemmatimonadetes Filterer 4 0.000000000 0.00000000 0.000000000
## 58 Gemmatimonadetes Predator 6 0.088888889 0.11863420 0.048432210
## 59 Gemmatimonadetes Scraper 9 0.000000000 0.00000000 0.000000000
## 60 Gemmatimonadetes Shredder 7 0.042857143 0.07382232 0.027902216
## 61 GN02 Filterer 4 0.633333333 1.13724814 0.568624070
## 62 GN02 Predator 6 0.083333333 0.11303883 0.046147910
## 63 GN02 Scraper 9 0.000000000 0.00000000 0.000000000
## 64 GN02 Shredder 7 0.009523810 0.01626500 0.006147593
## 65 Nitrospirae Filterer 4 0.000000000 0.00000000 0.000000000
## 66 Nitrospirae Predator 6 0.055555556 0.10680547 0.043603149
## 67 Nitrospirae Scraper 9 0.000000000 0.00000000 0.000000000
## 68 Nitrospirae Shredder 7 0.000000000 0.00000000 0.000000000
## 69 OD1 Filterer 4 0.016666667 0.03333333 0.016666667
## 70 OD1 Predator 6 0.022222222 0.04036867 0.016480441
## 71 OD1 Scraper 9 0.014814815 0.04444444 0.014814815
## 72 OD1 Shredder 7 0.004761905 0.01259882 0.004761905
## 73 Planctomycetes Filterer 4 1.308333333 0.40858836 0.204294178
## 74 Planctomycetes Predator 6 1.600000000 0.76361712 0.311745385
## 75 Planctomycetes Scraper 9 2.444444444 2.76551181 0.921837269
## 76 Planctomycetes Shredder 7 0.290476190 0.20522281 0.077566932
## 77 Proteobacteria Filterer 4 45.125000000 15.08895537 7.544477684
## 78 Proteobacteria Predator 6 51.183333333 9.36040716 3.821370222
## 79 Proteobacteria Scraper 9 68.003703704 27.10068903 9.033563009
## 80 Proteobacteria Shredder 7 35.666666667 9.46834883 3.578699477
## 81 SR1 Filterer 4 0.408333333 0.59961407 0.299807037
## 82 SR1 Predator 6 0.144444444 0.14246507 0.058161121
## 83 SR1 Scraper 9 0.003703704 0.01111111 0.003703704
## 84 SR1 Shredder 7 0.300000000 0.31856158 0.120404960
## 85 Synergistetes Filterer 4 0.600000000 0.89566859 0.447834295
## 86 Synergistetes Predator 6 0.000000000 0.00000000 0.000000000
## 87 Synergistetes Scraper 9 0.000000000 0.00000000 0.000000000
## 88 Synergistetes Shredder 7 0.080952381 0.13313477 0.050320214
## 89 Tenericutes Filterer 4 3.075000000 3.13479576 1.567397879
## 90 Tenericutes Predator 6 5.066666667 10.49615803 4.285038571
## 91 Tenericutes Scraper 9 24.988888889 28.21126646 9.403755488
## 92 Tenericutes Shredder 7 0.490476190 0.32071349 0.121218305
## 93 TM6 Filterer 4 0.000000000 0.00000000 0.000000000
## 94 TM6 Predator 6 0.011111111 0.02721655 0.011111111
## 95 TM6 Scraper 9 0.018518519 0.05555556 0.018518519
## 96 TM6 Shredder 7 0.000000000 0.00000000 0.000000000
## 97 TM7 Filterer 4 0.025000000 0.05000000 0.025000000
## 98 TM7 Predator 6 0.005555556 0.01360828 0.005555556
## 99 TM7 Scraper 9 0.081481481 0.23220627 0.077402091
## 100 TM7 Shredder 7 0.004761905 0.01259882 0.004761905
## 101 Verrucomicrobia Filterer 4 1.783333333 1.38684294 0.693421469
## 102 Verrucomicrobia Predator 6 3.827777778 2.46949538 1.008167265
## 103 Verrucomicrobia Scraper 9 0.740740741 0.56119361 0.187064538
## 104 Verrucomicrobia Shredder 7 0.409523810 0.16410733 0.062026742cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)
cdataplot#+ scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))
#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(
Trtdata2 <- ddply(df, c("Phylum", "FFG","Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
dfStats<-subset(df, FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station,scales = "free_x")
cdataplot2
########
compare_means(Abundance ~ FFG, data = df, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.0147 0.021 0.01466 * Kruskal-Wall~
## 2 Tenericutes Abundance 0.0643 0.075 0.06432 ns Kruskal-Wall~
## 3 Bacteroidetes Abundance 0.000400 0.0014 0.00040 *** Kruskal-Wall~
## 4 Firmicutes Abundance 0.000209 0.0014 0.00021 *** Kruskal-Wall~
## 5 Planctomycetes Abundance 0.0127 0.021 0.01274 * Kruskal-Wall~
## 6 Verrucomicrobia Abundance 0.00184 0.0043 0.00184 ** Kruskal-Wall~
## 7 Actinobacteria Abundance 0.118 0.12 0.11762 ns Kruskal-Wall~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Phylum", p.adjust.method = "fdr")
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Phylum))
for (i in levels(Means$Phylum)){
Tax<-i
TaxAbundance<-subset(Means,Phylum==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+10,label=vec)) +scale_fill_manual(values=cbPalette)
cdataplot
head(Trtdata)## Phylum FFG N mean sd se
## 1 Actinobacteria Filterer 4 0.3833333 0.4811252 0.2405626
## 2 Actinobacteria Predator 6 1.3111111 0.3862450 0.1576839
## 3 Actinobacteria Scraper 9 1.2259259 1.2891403 0.4297134
## 4 Actinobacteria Shredder 7 1.2476190 1.4698018 0.5555329
## 5 Bacteroidetes Filterer 4 22.6166667 3.7669616 1.8834808
## 6 Bacteroidetes Predator 6 32.7000000 10.7146006 4.3742174dfStatsPlot<-subset(Trtdata,Phylum!="Actinobacteria"&Phylum!= "Verrucomicrobia")
#vec<-c("a","b","a","b","a","b","a","b","a","b")
dfStatsPlot # Only show significant taxa ## Phylum FFG N mean sd se
## 5 Bacteroidetes Filterer 4 22.6166667 3.7669616 1.88348082
## 6 Bacteroidetes Predator 6 32.7000000 10.7146006 4.37421739
## 7 Bacteroidetes Scraper 9 0.8444444 0.6502136 0.21673788
## 8 Bacteroidetes Shredder 7 22.8380952 6.0532033 2.28789578
## 9 Firmicutes Filterer 4 20.0250000 13.0079285 6.50396425
## 10 Firmicutes Predator 6 0.2555556 0.4385160 0.17902341
## 11 Firmicutes Scraper 9 0.7888889 1.2328828 0.41096093
## 12 Firmicutes Shredder 7 36.6333333 7.7272486 2.92062543
## 13 Planctomycetes Filterer 4 1.3083333 0.4085884 0.20429418
## 14 Planctomycetes Predator 6 1.6000000 0.7636171 0.31174539
## 15 Planctomycetes Scraper 9 2.4444444 2.7655118 0.92183727
## 16 Planctomycetes Shredder 7 0.2904762 0.2052228 0.07756693
## 17 Proteobacteria Filterer 4 45.1250000 15.0889554 7.54447768
## 18 Proteobacteria Predator 6 51.1833333 9.3604072 3.82137022
## 19 Proteobacteria Scraper 9 68.0037037 27.1006890 9.03356301
## 20 Proteobacteria Shredder 7 35.6666667 9.4683488 3.57869948
## 21 Tenericutes Filterer 4 3.0750000 3.1347958 1.56739788
## 22 Tenericutes Predator 6 5.0666667 10.4961580 4.28503857
## 23 Tenericutes Scraper 9 24.9888889 28.2112665 9.40375549
## 24 Tenericutes Shredder 7 0.4904762 0.3207135 0.12121831cdataplot=ggplot(dfStatsPlot, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+#geom_text(aes(x=FFG, y=mean+se+10,label=vec))+
scale_fill_manual(values=cbPalette)+theme(legend.justification=c(0.05,0.95), legend.position=c(0.7,0.4))
cdataplot
#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Phylum))
#SigList<-length(unique(Trtdata$Phylum))
#SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")
Means## # A tibble: 7 x 9
## Phylum .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobac~ Abunda~ Scraper Shred~ 0.0149 0.021 0.015 * Wilco~
## 2 Tenericut~ Abunda~ Scraper Shred~ 0.0148 0.021 0.015 * Wilco~
## 3 Firmicutes Abunda~ Scraper Shred~ 0.00103 0.0036 0.001 ** Wilco~
## 4 Bacteroid~ Abunda~ Scraper Shred~ 0.00103 0.0036 0.001 ** Wilco~
## 5 Planctomy~ Abunda~ Scraper Shred~ 0.0148 0.021 0.015 * Wilco~
## 6 Actinobac~ Abunda~ Scraper Shred~ 0.751 0.75 0.751 ns Wilco~
## 7 Verrucomi~ Abunda~ Scraper Shred~ 0.289 0.34 0.289 ns Wilco~MeansStation## # A tibble: 7 x 9
## Phylum .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacte~ Abunda~ Ostana PDR 0.958 1 0.958 ns Wilcox~
## 2 Tenericutes Abunda~ Ostana PDR 0.528 0.89 0.528 ns Wilcox~
## 3 Firmicutes Abunda~ Ostana PDR 1 1 1.000 ns Wilcox~
## 4 Bacteroidet~ Abunda~ Ostana PDR 0.637 0.89 0.637 ns Wilcox~
## 5 Planctomyce~ Abunda~ Ostana PDR 0.0180 0.13 0.018 * Wilcox~
## 6 Actinobacte~ Abunda~ Ostana PDR 0.495 0.89 0.495 ns Wilcox~
## 7 Verrucomicr~ Abunda~ Ostana PDR 0.0738 0.26 0.074 ns Wilcox~Means<-compare_means(Abundance ~ FFG, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")
Trtdata <- ddply(dfStats, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Phylum))
for (i in levels(Means$Phylum)){
Tax<-i
TaxAbundance<-subset(Means,Phylum==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+5,label=vec))+scale_fill_manual(values=cbPalette)
cdataplot
#######Split by Station
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")
dfStatsOstana<-subset(dfOstana,FFG!="Predator")
print("Ostana")## [1] "Ostana"compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.291 0.290 0.2906 ns Kruskal-Wallis
## 2 Tenericutes Abundance 0.0549 0.077 0.0549 ns Kruskal-Wallis
## 3 Firmicutes Abundance 0.0121 0.042 0.0121 * Kruskal-Wallis
## 4 Bacteroidetes Abundance 0.00990 0.042 0.0099 ** Kruskal-Wallis
## 5 Planctomycetes Abundance 0.0430 0.077 0.0430 * Kruskal-Wallis
## 6 Actinobacteria Abundance 0.109 0.13 0.1089 ns Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.0484 0.077 0.0484 * Kruskal-WallisMeans<-compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## [1] Phylum group1 group2 p.adj
## <0 rows> (or 0-length row.names)print("Pian della Regina")## [1] "Pian della Regina"Means<-compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Phylum group1 group2 p.adj
## 2 Proteobacteria Predator Shredder 0.071
## 3 Proteobacteria Scraper Shredder 0.071
## 7 Bacteroidetes Predator Scraper 0.071
## 9 Bacteroidetes Scraper Shredder 0.071
## 11 Firmicutes Predator Shredder 0.071
## 12 Firmicutes Scraper Shredder 0.071
## 13 Verrucomicrobia Predator Scraper 0.071
## 14 Verrucomicrobia Predator Shredder 0.071
## 20 Planctomycetes Predator Shredder 0.071#########
############
ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Ostana(>1%)")
ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance PDR(>1%)")
print("Family level")## [1] "Family level"df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")
dfStatsOstana<-subset(dfOstana,FFG!="Predator")
print("Ostana")## [1] "Ostana"compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aeromonadaceae Abundance 0.0817 0.091 0.0817 ns Kruskal-Wal~
## 2 Comamonadaceae Abundance 0.00786 0.021 0.0079 ** Kruskal-Wal~
## 3 Rhodobacteraceae Abundance 0.886 0.89 0.8865 ns Kruskal-Wal~
## 4 Ruminococcaceae Abundance 0.00712 0.021 0.0071 ** Kruskal-Wal~
## 5 Desulfovibrionac~ Abundance 0.0105 0.021 0.0105 * Kruskal-Wal~
## 6 Lachnospiraceae Abundance 0.0118 0.021 0.0118 * Kruskal-Wal~
## 7 Rikenellaceae Abundance 0.0330 0.041 0.0330 * Kruskal-Wal~
## 8 Methylophilaceae Abundance 0.0139 0.021 0.0139 * Kruskal-Wal~
## 9 Enterobacteriace~ Abundance 0.0144 0.021 0.0144 * Kruskal-Wal~
## 10 Cytophagaceae Abundance 0.0139 0.021 0.0139 * Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Family group1 group2 p.adj
## 3 Comamonadaceae Filterer Scraper 0.058
## 4 Comamonadaceae Filterer Shredder 0.088
## 5 Comamonadaceae Scraper Shredder 0.076
## 9 Ruminococcaceae Filterer Scraper 0.054
## 11 Ruminococcaceae Scraper Shredder 0.054
## 12 Desulfovibrionaceae Filterer Scraper 0.054
## 14 Desulfovibrionaceae Scraper Shredder 0.054
## 15 Lachnospiraceae Filterer Scraper 0.054
## 17 Lachnospiraceae Scraper Shredder 0.054
## 18 Rikenellaceae Filterer Scraper 0.079
## 20 Rikenellaceae Scraper Shredder 0.054
## 21 Methylophilaceae Filterer Scraper 0.054
## 23 Methylophilaceae Scraper Shredder 0.063
## 25 Enterobacteriaceae Filterer Shredder 0.079
## 26 Enterobacteriaceae Scraper Shredder 0.076
## 27 Cytophagaceae Filterer Scraper 0.054
## 29 Cytophagaceae Scraper Shredder 0.063print("Pian della Regina")## [1] "Pian della Regina"compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteriace~ Abundance 0.0366 0.041 0.0366 * Kruskal-Wal~
## 2 Cytophagaceae Abundance 0.00609 0.013 0.0061 ** Kruskal-Wal~
## 3 Desulfovibrionac~ Abundance 0.00537 0.013 0.0054 ** Kruskal-Wal~
## 4 Lachnospiraceae Abundance 0.00921 0.013 0.0092 ** Kruskal-Wal~
## 5 Rikenellaceae Abundance 0.00537 0.013 0.0054 ** Kruskal-Wal~
## 6 Rhodobacteraceae Abundance 0.0249 0.031 0.0249 * Kruskal-Wal~
## 7 Comamonadaceae Abundance 0.00715 0.013 0.0072 ** Kruskal-Wal~
## 8 Ruminococcaceae Abundance 0.00921 0.013 0.0092 ** Kruskal-Wal~
## 9 Aeromonadaceae Abundance 0.368 0.37 0.3679 ns Kruskal-Wal~
## 10 Methylophilaceae Abundance 0.00609 0.013 0.0061 ** Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Family group1 group2 p.adj
## 3 Enterobacteriaceae Scraper Shredder 0.041
## 4 Cytophagaceae Predator Scraper 0.041
## 5 Cytophagaceae Predator Shredder 0.041
## 6 Cytophagaceae Scraper Shredder 0.041
## 7 Desulfovibrionaceae Predator Shredder 0.041
## 8 Desulfovibrionaceae Scraper Shredder 0.041
## 10 Lachnospiraceae Predator Shredder 0.041
## 11 Lachnospiraceae Scraper Shredder 0.041
## 12 Rikenellaceae Predator Shredder 0.041
## 13 Rikenellaceae Scraper Shredder 0.041
## 14 Rhodobacteraceae Predator Scraper 0.041
## 15 Rhodobacteraceae Predator Shredder 0.041
## 17 Comamonadaceae Predator Scraper 0.041
## 18 Comamonadaceae Predator Shredder 0.041
## 19 Comamonadaceae Scraper Shredder 0.041
## 21 Ruminococcaceae Predator Shredder 0.041
## 22 Ruminococcaceae Scraper Shredder 0.041
## 25 Methylophilaceae Predator Scraper 0.041
## 26 Methylophilaceae Predator Shredder 0.041
## 27 Methylophilaceae Scraper Shredder 0.041ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Ostana(>2%)")
ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance PDR(>2%)")
print("Genus Level")## [1] "Genus Level"df <- psmelt(GenusLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")
dfStatsOstana<-subset(dfOstana,FFG!="Predator")
print("Ostana")## [1] "Ostana"compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 3 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacter Abundance 0.886 0.89 0.886 ns Kruskal-Wallis
## 2 Clostridium Abundance 0.0118 0.021 0.012 * Kruskal-Wallis
## 3 Methylotenera Abundance 0.0139 0.021 0.014 * Kruskal-Walliscompare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Genus", p.adjust.method = "fdr")## # A tibble: 9 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobact~ Abunda~ Filter~ Scraper 0.903 0.9 0.903 ns Wilcox~
## 2 Rhodobact~ Abunda~ Filter~ Shredd~ 0.860 0.9 0.860 ns Wilcox~
## 3 Rhodobact~ Abunda~ Scraper Shredd~ 0.766 0.9 0.766 ns Wilcox~
## 4 Clostridi~ Abunda~ Filter~ Scraper 0.0108 0.05 0.011 * Wilcox~
## 5 Clostridi~ Abunda~ Filter~ Shredd~ 0.596 0.89 0.596 ns Wilcox~
## 6 Clostridi~ Abunda~ Scraper Shredd~ 0.0168 0.05 0.017 * Wilcox~
## 7 Methylote~ Abunda~ Filter~ Scraper 0.0151 0.05 0.015 * Wilcox~
## 8 Methylote~ Abunda~ Filter~ Shredd~ 0.596 0.89 0.596 ns Wilcox~
## 9 Methylote~ Abunda~ Scraper Shredd~ 0.0262 0.059 0.026 * Wilcox~print("Pian della Regina")## [1] "Pian della Regina"compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 3 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Clostridium Abundance 0.00537 0.0091 0.0054 ** Kruskal-Wallis
## 2 Rhodobacter Abundance 0.0246 0.025 0.0246 * Kruskal-Wallis
## 3 Methylotenera Abundance 0.00609 0.0091 0.0061 ** Kruskal-Walliscompare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Genus", p.adjust.method = "fdr")## # A tibble: 8 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Clostridi~ Abunda~ Predat~ Shredd~ 0.0211 0.035 0.021 * Wilcox~
## 2 Clostridi~ Abunda~ Scraper Shredd~ 0.0211 0.035 0.021 * Wilcox~
## 3 Rhodobact~ Abunda~ Predat~ Scraper 0.0304 0.035 0.030 * Wilcox~
## 4 Rhodobact~ Abunda~ Predat~ Shredd~ 0.0294 0.035 0.029 * Wilcox~
## 5 Rhodobact~ Abunda~ Scraper Shredd~ 1 1 1.000 ns Wilcox~
## 6 Methylote~ Abunda~ Predat~ Scraper 0.0211 0.035 0.021 * Wilcox~
## 7 Methylote~ Abunda~ Predat~ Shredd~ 0.0304 0.035 0.030 * Wilcox~
## 8 Methylote~ Abunda~ Scraper Shredd~ 0.0211 0.035 0.021 * Wilcox~ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Genus)+ylab("Relative Abundance Ostana(>2%)")
ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Genus)+ylab("Relative Abundance PDR(>2%)")
Split by FFG
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
dfShredders<-subset(df,FFG =="Shredder")
dfScrapers<-subset(df,FFG == "Scraper")
print("Shredders")## [1] "Shredders"compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.480 0.56 0.480 ns Kruskal-Wallis
## 2 Firmicutes Abundance 0.480 0.56 0.480 ns Kruskal-Wallis
## 3 Bacteroidetes Abundance 0.0339 0.18 0.034 * Kruskal-Wallis
## 4 Actinobacteria Abundance 0.0771 0.18 0.077 ns Kruskal-Wallis
## 5 Tenericutes Abundance 0.372 0.56 0.372 ns Kruskal-Wallis
## 6 Planctomycetes Abundance 0.0771 0.18 0.077 ns Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.589 0.59 0.589 ns Kruskal-WallisMeans<-compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Phylum group1 group2 p
## 3 Bacteroidetes Ostana PDR 0.05182993print("Scrapers")## [1] "Scrapers"compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.327 0.570 0.33 ns Kruskal-Wallis
## 2 Tenericutes Abundance 0.462 0.62 0.46 ns Kruskal-Wallis
## 3 Planctomycetes Abundance 0.0500 0.35 0.05 ns Kruskal-Wallis
## 4 Actinobacteria Abundance 0.623 0.62 0.62 ns Kruskal-Wallis
## 5 Firmicutes Abundance 0.623 0.62 0.62 ns Kruskal-Wallis
## 6 Bacteroidetes Abundance 0.221 0.51 0.22 ns Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.142 0.5 0.14 ns Kruskal-WallisMeans<-compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Phylum", p.adjust.method = "fdr")
keeps <- c("Phylum","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Phylum group1 group2 p
## 3 Planctomycetes Ostana PDR 0.06619258ggplot(dfShredders, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Shredders(>1%)")
ggplot(dfScrapers, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Scrapers (>1%)")
############
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
dfShredders<-subset(df,FFG =="Shredder")
dfScrapers<-subset(df,FFG == "Scraper")
Means<-compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Family group1 group2 p
## 3 Lachnospiraceae Ostana PDR 0.05182993
## 4 Rikenellaceae Ostana PDR 0.05182993
## 6 Rhodobacteraceae Ostana PDR 0.05182993
## 7 Enterobacteriaceae Ostana PDR 0.05182993print("Scrapers")## [1] "Scrapers"compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteriac~ Abundan~ 0.0143 0.086 0.014 * Kruskal-Wa~
## 2 Aeromonadaceae Abundan~ 0.283 0.37 0.283 ns Kruskal-Wa~
## 3 Rhodobacteraceae Abundan~ 0.142 0.34 0.142 ns Kruskal-Wa~
## 4 Comamonadaceae Abundan~ 0.169 0.34 0.169 ns Kruskal-Wa~
## 5 Methylophilaceae Abundan~ 0.371 0.37 0.371 ns Kruskal-Wa~
## 6 Cytophagaceae Abundan~ 0.371 0.37 0.371 ns Kruskal-Wa~
## 7 Desulfovibriona~ Abundan~ NaN NaN NA ? Kruskal-Wa~
## 8 Rikenellaceae Abundan~ NaN NaN NA ? Kruskal-Wa~
## 9 Ruminococcaceae Abundan~ NaN NaN NA ? Kruskal-Wa~
## 10 Lachnospiraceae Abundan~ NaN NaN NA ? Kruskal-Wa~Means<-compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p","method","p.signif")
keeps=Means[keeps]
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]
FilteredResults## Family group1 group2 p
## 1 Enterobacteriaceae Ostana PDR 0.01996445ggplot(dfShredders, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Shredders(>1%)")
ggplot(dfScrapers, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Scrapers (>1%)")
Family Level
The top 5 Families across all samples are given in the table below
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
dfStats<-subset(df,FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
df2 <- psmelt(FamilyAll)
df2$Abundance<-df2$Abundance*100
FamilyAll<-ddply(df2, c("Family"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
FamilySorted<-FamilyAll[order(-FamilyAll$mean),]
#FamilySorted[1:5,]
FamilyAll2<-ddply(df2, c("Family","FFG"), summarise, #For Rel abu tabs
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
FamilySorted2<-FamilyAll2[order(-FamilyAll$mean),]
FamilySorted2## Family FFG N mean sd se
## 62 Aeromonadaceae Predator 6 0.000000000 0.00000000 0.000000000
## 56 Acidobacteriaceae Shredder 7 0.004761905 0.01259882 0.004761905
## 126 Brevibacteriaceae Predator 6 0.038888889 0.09525793 0.038888889
## 43 A1-B1 Scraper 9 0.000000000 0.00000000 0.000000000
## 49 Acetobacteraceae Filterer 4 0.050000000 0.05773503 0.028867513
## 133 Burkholderiaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 84 Anaeroplasmataceae Shredder 7 0.090476190 0.08758971 0.033105799
## 16 [Exiguobacteraceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 132 Brucellaceae Shredder 7 0.014285714 0.02622653 0.009912695
## 90 Armatimonadaceae Predator 6 0.155555556 0.15587269 0.063634760
## 135 Burkholderiaceae Scraper 9 0.029629630 0.08888889 0.029629630
## 119 Bifidobacteriaceae Scraper 9 0.003703704 0.01111111 0.003703704
## 149 Chamaesiphonaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 95 auto67_4W Scraper 9 0.000000000 0.00000000 0.000000000
## 33 0319-6G20 Filterer 4 0.000000000 0.00000000 0.000000000
## 39 211ds20 Scraper 9 0.000000000 0.00000000 0.000000000
## 69 Alcaligenaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 147 Caulobacteraceae Scraper 9 0.125925926 0.37777778 0.125925926
## 89 Armatimonadaceae Filterer 4 0.016666667 0.01924501 0.009622504
## 117 Bifidobacteriaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 50 Acetobacteraceae Predator 6 0.016666667 0.02788867 0.011385501
## 114 Beijerinckiaceae Predator 6 0.000000000 0.00000000 0.000000000
## 150 Chamaesiphonaceae Predator 6 0.027777778 0.05340273 0.021801574
## 141 Caldicoprobacteraceae Filterer 4 0.000000000 0.00000000 0.000000000
## 98 Bacillaceae Predator 6 0.011111111 0.02721655 0.011111111
## 78 Anaerobrancaceae Predator 6 0.000000000 0.00000000 0.000000000
## 115 Beijerinckiaceae Scraper 9 0.400000000 0.80484643 0.268282144
## 146 Caulobacteraceae Predator 6 0.277777778 0.13109228 0.053518198
## 113 Beijerinckiaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 128 Brevibacteriaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 54 Acidobacteriaceae Predator 6 0.000000000 0.00000000 0.000000000
## 123 Bradyrhizobiaceae Scraper 9 0.122222222 0.25927249 0.086424162
## 73 Alteromonadaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 129 Brucellaceae Filterer 4 0.008333333 0.01666667 0.008333333
## 86 Anaplasmataceae Predator 6 0.027777778 0.06804138 0.027777778
## 3 [Bryobacteraceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 104 Bacteriovoracaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 64 Aeromonadaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 29 [Weeksellaceae] Filterer 4 0.008333333 0.01666667 0.008333333
## 48 A4b Shredder 7 0.000000000 0.00000000 0.000000000
## 37 211ds20 Filterer 4 0.000000000 0.00000000 0.000000000
## 41 A1-B1 Filterer 4 0.000000000 0.00000000 0.000000000
## 28 [Tissierellaceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 42 A1-B1 Predator 6 0.038888889 0.06116160 0.024969117
## 106 Bacteroidaceae Predator 6 0.000000000 0.00000000 0.000000000
## 125 Brevibacteriaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 61 Aeromonadaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 75 Alteromonadaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 6 [Cerasicoccaceae] Predator 6 0.005555556 0.01360828 0.005555556
## 8 [Cerasicoccaceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 140 C111 Shredder 7 0.019047619 0.01781742 0.006734350
## 25 [Tissierellaceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 107 Bacteroidaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 143 Caldicoprobacteraceae Scraper 9 0.000000000 0.00000000 0.000000000
## 118 Bifidobacteriaceae Predator 6 0.000000000 0.00000000 0.000000000
## 57 Actinomycetaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 124 Bradyrhizobiaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 138 C111 Predator 6 0.150000000 0.16295875 0.066527633
## 68 AK1AB1_02E Shredder 7 0.000000000 0.00000000 0.000000000
## 92 Armatimonadaceae Shredder 7 0.095238095 0.11126973 0.042056004
## 26 [Tissierellaceae] Predator 6 0.000000000 0.00000000 0.000000000
## 35 0319-6G20 Scraper 9 0.000000000 0.00000000 0.000000000
## 23 [Mogibacteriaceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 102 Bacteriovoracaceae Predator 6 0.283333333 0.14414499 0.058846945
## 32 [Weeksellaceae] Shredder 7 0.085714286 0.09786072 0.036987874
## 45 A4b Filterer 4 0.000000000 0.00000000 0.000000000
## 139 C111 Scraper 9 0.040740741 0.12222222 0.040740741
## 13 [Exiguobacteraceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 20 [Fimbriimonadaceae] Shredder 7 0.009523810 0.01626500 0.006147593
## 31 [Weeksellaceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 24 [Mogibacteriaceae] Shredder 7 0.352380952 0.36201961 0.136830553
## 44 A1-B1 Shredder 7 0.042857143 0.07382232 0.027902216
## 58 Actinomycetaceae Predator 6 0.000000000 0.00000000 0.000000000
## 101 Bacteriovoracaceae Filterer 4 0.016666667 0.03333333 0.016666667
## 105 Bacteroidaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 144 Caldicoprobacteraceae Shredder 7 0.028571429 0.04049953 0.015307382
## 19 [Fimbriimonadaceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 91 Armatimonadaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 110 Bdellovibrionaceae Predator 6 0.244444444 0.22377237 0.091354688
## 122 Bradyrhizobiaceae Predator 6 0.005555556 0.01360828 0.005555556
## 53 Acidobacteriaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 38 211ds20 Predator 6 0.005555556 0.01360828 0.005555556
## 2 [Bryobacteraceae] Predator 6 0.083333333 0.11105554 0.045338235
## 5 [Cerasicoccaceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 134 Burkholderiaceae Predator 6 0.000000000 0.00000000 0.000000000
## 77 Anaerobrancaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 21 [Mogibacteriaceae] Filterer 4 0.025000000 0.03191424 0.015957118
## 52 Acetobacteraceae Shredder 7 0.023809524 0.04178554 0.015793451
## 103 Bacteriovoracaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 99 Bacillaceae Scraper 9 0.237037037 0.71111111 0.237037037
## 81 Anaeroplasmataceae Filterer 4 0.000000000 0.00000000 0.000000000
## 142 Caldicoprobacteraceae Predator 6 0.000000000 0.00000000 0.000000000
## 11 [Chthoniobacteraceae] Scraper 9 0.166666667 0.28037673 0.093458910
## 130 Brucellaceae Predator 6 0.066666667 0.06666667 0.027216553
## 1 [Bryobacteraceae] Filterer 4 0.000000000 0.00000000 0.000000000
## 72 Alcaligenaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 120 Bifidobacteriaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 46 A4b Predator 6 0.022222222 0.04036867 0.016480441
## 63 Aeromonadaceae Scraper 9 8.733333333 15.74738955 5.249129851
## 74 Alteromonadaceae Predator 6 0.111111111 0.22673202 0.092562958
## 97 Bacillaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 14 [Exiguobacteraceae] Predator 6 0.005555556 0.01360828 0.005555556
## 85 Anaplasmataceae Filterer 4 0.000000000 0.00000000 0.000000000
## 112 Bdellovibrionaceae Shredder 7 0.028571429 0.05245305 0.019825390
## 79 Anaerobrancaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 100 Bacillaceae Shredder 7 0.009523810 0.02519763 0.009523810
## 18 [Fimbriimonadaceae] Predator 6 0.100000000 0.11737878 0.047919686
## 34 0319-6G20 Predator 6 0.016666667 0.04082483 0.016666667
## 65 AK1AB1_02E Filterer 4 0.008333333 0.01666667 0.008333333
## 36 0319-6G20 Shredder 7 0.000000000 0.00000000 0.000000000
## 121 Bradyrhizobiaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 22 [Mogibacteriaceae] Predator 6 0.000000000 0.00000000 0.000000000
## 83 Anaeroplasmataceae Scraper 9 0.000000000 0.00000000 0.000000000
## 82 Anaeroplasmataceae Predator 6 0.000000000 0.00000000 0.000000000
## 94 auto67_4W Predator 6 0.150000000 0.16158933 0.065968567
## 131 Brucellaceae Scraper 9 4.648148148 6.95218014 2.317393379
## 12 [Chthoniobacteraceae] Shredder 7 0.038095238 0.06214848 0.023489918
## 59 Actinomycetaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 109 Bdellovibrionaceae Filterer 4 0.158333333 0.14240006 0.071200031
## 9 [Chthoniobacteraceae] Filterer 4 0.275000000 0.18534253 0.092671263
## 66 AK1AB1_02E Predator 6 0.000000000 0.00000000 0.000000000
## 67 AK1AB1_02E Scraper 9 0.000000000 0.00000000 0.000000000
## 127 Brevibacteriaceae Scraper 9 0.129629630 0.38888889 0.129629630
## 148 Caulobacteraceae Shredder 7 0.085714286 0.22677868 0.085714286
## 88 Anaplasmataceae Shredder 7 0.000000000 0.00000000 0.000000000
## 145 Caulobacteraceae Filterer 4 0.016666667 0.03333333 0.016666667
## 15 [Exiguobacteraceae] Scraper 9 0.000000000 0.00000000 0.000000000
## 40 211ds20 Shredder 7 0.000000000 0.00000000 0.000000000
## 55 Acidobacteriaceae Scraper 9 0.037037037 0.11111111 0.037037037
## 60 Actinomycetaceae Shredder 7 0.009523810 0.02519763 0.009523810
## 87 Anaplasmataceae Scraper 9 0.000000000 0.00000000 0.000000000
## 93 auto67_4W Filterer 4 0.050000000 0.05773503 0.028867513
## 96 auto67_4W Shredder 7 0.000000000 0.00000000 0.000000000
## 108 Bacteroidaceae Shredder 7 0.004761905 0.01259882 0.004761905
## 111 Bdellovibrionaceae Scraper 9 0.114814815 0.23339946 0.077799821
## 136 Burkholderiaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 137 C111 Filterer 4 0.091666667 0.18333333 0.091666667
## 4 [Bryobacteraceae] Shredder 7 0.000000000 0.00000000 0.000000000
## 7 [Cerasicoccaceae] Scraper 9 0.081481481 0.18189673 0.060632243
## 10 [Chthoniobacteraceae] Predator 6 0.333333333 0.36757463 0.150061716
## 17 [Fimbriimonadaceae] Filterer 4 0.016666667 0.03333333 0.016666667
## 27 [Tissierellaceae] Scraper 9 0.003703704 0.01111111 0.003703704
## 30 [Weeksellaceae] Predator 6 0.300000000 0.35402448 0.144529889
## 47 A4b Scraper 9 0.000000000 0.00000000 0.000000000
## 51 Acetobacteraceae Scraper 9 0.088888889 0.16914819 0.056382731
## 70 Alcaligenaceae Predator 6 0.000000000 0.00000000 0.000000000
## 71 Alcaligenaceae Scraper 9 0.029629630 0.07718024 0.025726748
## 76 Alteromonadaceae Shredder 7 0.042857143 0.11338934 0.042857143
## 80 Anaerobrancaceae Shredder 7 0.161904762 0.30088054 0.113722155
## 116 Beijerinckiaceae Shredder 7 0.000000000 0.00000000 0.000000000cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 2%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)
cdataplot#+ scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))
#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(
Trtdata2 <- ddply(df, c("Family", "FFG","Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station)
cdataplot2
########
dfStats<-subset(df, FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
compare_means(Abundance ~ FFG, data = dfStats, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")# Differences between shredder and scrapers only## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteriac~ Abundan~ 0.340 3.80e-1 0.34041 ns Kruskal-Wa~
## 2 Aeromonadaceae Abundan~ 0.0516 6.40e-2 0.05155 ns Kruskal-Wa~
## 3 Rhodobacteraceae Abundan~ 0.958 9.60e-1 0.95779 ns Kruskal-Wa~
## 4 Ruminococcaceae Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 5 Desulfovibriona~ Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 6 Lachnospiraceae Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 7 Rikenellaceae Abundan~ 0.000239 6.00e-4 0.00024 *** Kruskal-Wa~
## 8 Methylophilaceae Abundan~ 0.000369 6.10e-4 0.00037 *** Kruskal-Wa~
## 9 Comamonadaceae Abundan~ 0.000818 1.20e-3 0.00082 *** Kruskal-Wa~
## 10 Cytophagaceae Abundan~ 0.000369 6.10e-4 0.00037 *** Kruskal-Wa~#df<-subset(df, Family!="Aeromonadaceae"&Family!="Rhodobacteraceae"&Family!="Enterobacteriaceae")# Remove groups that were not sig by Kruksal-Wallis
#df <- df[which(df$Family!="Aeromonadaceae")] #& df$Family!="Rhodobacteraceae"& df$Family!="Enterobacteriaceae"
#df
Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Family", p.adjust.method = "fdr")
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Family))
for (i in levels(Means$Family)){
Tax<-i
TaxAbundance<-subset(Means,Family==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+5,label=vec))+scale_fill_manual(values=cbPalette)
cdataplot
(Trtdata)## Family FFG N mean sd se
## 1 Aeromonadaceae Filterer 4 0.000000000 0.00000000 0.000000000
## 2 Aeromonadaceae Predator 6 0.000000000 0.00000000 0.000000000
## 3 Aeromonadaceae Scraper 9 8.733333333 15.74738955 5.249129851
## 4 Aeromonadaceae Shredder 7 0.000000000 0.00000000 0.000000000
## 5 Comamonadaceae Filterer 4 12.333333333 8.80151502 4.400757511
## 6 Comamonadaceae Predator 6 6.294444444 2.43715833 0.994965723
## 7 Comamonadaceae Scraper 9 0.362962963 0.46050831 0.153502769
## 8 Comamonadaceae Shredder 7 2.314285714 1.02068553 0.385782867
## 9 Cytophagaceae Filterer 4 1.500000000 1.09273696 0.546368482
## 10 Cytophagaceae Predator 6 14.833333333 8.29572849 3.386716973
## 11 Cytophagaceae Scraper 9 0.011111111 0.03333333 0.011111111
## 12 Cytophagaceae Shredder 7 1.171428571 0.61626971 0.232928057
## 13 Desulfovibrionaceae Filterer 4 10.708333333 5.08821259 2.544106297
## 14 Desulfovibrionaceae Predator 6 0.255555556 0.62598071 0.255555556
## 15 Desulfovibrionaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 16 Desulfovibrionaceae Shredder 7 9.923809524 3.25671468 1.230922446
## 17 Enterobacteriaceae Filterer 4 0.025000000 0.05000000 0.025000000
## 18 Enterobacteriaceae Predator 6 8.122222222 15.38135475 6.279411783
## 19 Enterobacteriaceae Scraper 9 24.911111111 39.25302182 13.084340608
## 20 Enterobacteriaceae Shredder 7 1.985714286 2.21926390 0.838802912
## 21 Lachnospiraceae Filterer 4 6.000000000 5.05253878 2.526269391
## 22 Lachnospiraceae Predator 6 0.011111111 0.02721655 0.011111111
## 23 Lachnospiraceae Scraper 9 0.000000000 0.00000000 0.000000000
## 24 Lachnospiraceae Shredder 7 7.909523810 3.33042730 1.258783201
## 25 Methylophilaceae Filterer 4 2.508333333 1.45178689 0.725893447
## 26 Methylophilaceae Predator 6 3.155555556 3.02571693 1.235243766
## 27 Methylophilaceae Scraper 9 0.011111111 0.03333333 0.011111111
## 28 Methylophilaceae Shredder 7 5.090476190 2.19145768 0.828293148
## 29 Rhodobacteraceae Filterer 4 3.750000000 3.46297881 1.731489404
## 30 Rhodobacteraceae Predator 6 7.327777778 3.02591890 1.235326218
## 31 Rhodobacteraceae Scraper 9 3.977777778 6.42209727 2.140699090
## 32 Rhodobacteraceae Shredder 7 2.647619048 2.44198289 0.922982775
## 33 Rikenellaceae Filterer 4 4.725000000 4.26287809 2.131439046
## 34 Rikenellaceae Predator 6 0.000000000 0.00000000 0.000000000
## 35 Rikenellaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 36 Rikenellaceae Shredder 7 7.076190476 3.38677790 1.280081725
## 37 Ruminococcaceae Filterer 4 5.250000000 3.50518135 1.752590675
## 38 Ruminococcaceae Predator 6 0.005555556 0.01360828 0.005555556
## 39 Ruminococcaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 40 Ruminococcaceae Shredder 7 9.690476190 4.86722654 1.839638714TrtdataPlot<-subset(Trtdata, Family!="Aeromonadaceae"&Family!="Rhodobacteraceae"&Family!="Enterobacteriaceae"&Family!="Cytophagaceae")
(TrtdataPlot)## Family FFG N mean sd se
## 5 Comamonadaceae Filterer 4 12.333333333 8.80151502 4.400757511
## 6 Comamonadaceae Predator 6 6.294444444 2.43715833 0.994965723
## 7 Comamonadaceae Scraper 9 0.362962963 0.46050831 0.153502769
## 8 Comamonadaceae Shredder 7 2.314285714 1.02068553 0.385782867
## 13 Desulfovibrionaceae Filterer 4 10.708333333 5.08821259 2.544106297
## 14 Desulfovibrionaceae Predator 6 0.255555556 0.62598071 0.255555556
## 15 Desulfovibrionaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 16 Desulfovibrionaceae Shredder 7 9.923809524 3.25671468 1.230922446
## 21 Lachnospiraceae Filterer 4 6.000000000 5.05253878 2.526269391
## 22 Lachnospiraceae Predator 6 0.011111111 0.02721655 0.011111111
## 23 Lachnospiraceae Scraper 9 0.000000000 0.00000000 0.000000000
## 24 Lachnospiraceae Shredder 7 7.909523810 3.33042730 1.258783201
## 25 Methylophilaceae Filterer 4 2.508333333 1.45178689 0.725893447
## 26 Methylophilaceae Predator 6 3.155555556 3.02571693 1.235243766
## 27 Methylophilaceae Scraper 9 0.011111111 0.03333333 0.011111111
## 28 Methylophilaceae Shredder 7 5.090476190 2.19145768 0.828293148
## 33 Rikenellaceae Filterer 4 4.725000000 4.26287809 2.131439046
## 34 Rikenellaceae Predator 6 0.000000000 0.00000000 0.000000000
## 35 Rikenellaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 36 Rikenellaceae Shredder 7 7.076190476 3.38677790 1.280081725
## 37 Ruminococcaceae Filterer 4 5.250000000 3.50518135 1.752590675
## 38 Ruminococcaceae Predator 6 0.005555556 0.01360828 0.005555556
## 39 Ruminococcaceae Scraper 9 0.000000000 0.00000000 0.000000000
## 40 Ruminococcaceae Shredder 7 9.690476190 4.86722654 1.839638714#vec<-c("a","b","a","b","a","b","a","b","a","b","a","b")
cdataplot=ggplot(TrtdataPlot, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)#theme(axis.text.x = element_text(angle = 45, hjust = 1))+
cdataplot
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Family", p.adjust.method = "fdr")
#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Family))
#SigList<-length(unique(Trtdata$Family))
#SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
#Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Family", p.adjust.method = "fdr")
#for (i in levels(Means$Family)){
# Tax<-i
# TaxAbundance<-subset(Means,Family==i )
# Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
# difference<-TaxAbundance$p.adj
# names(difference)<-Hyphenated
# Letters<-multcompLetters(difference)
#print(Letters)
# SigList[i]<-Letters
#}
#vec<-unlist(SigList)
#vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
#+geom_text(aes(x=FFG, y=mean+se+5,label=vec))
Means=compare_means(Abundance ~ FFG, data = dfStats,
group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
#keeps
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]
FilteredResults## Family group1 group2 p.adj
## 2 Aeromonadaceae Scraper Shredder 0.07600
## 4 Ruminococcaceae Scraper Shredder 0.00075
## 5 Desulfovibrionaceae Scraper Shredder 0.00075
## 6 Lachnospiraceae Scraper Shredder 0.00075
## 7 Rikenellaceae Scraper Shredder 0.00075
## 8 Methylophilaceae Scraper Shredder 0.00076
## 9 Comamonadaceae Scraper Shredder 0.00140
## 10 Cytophagaceae Scraper Shredder 0.00076Genus level
The top 5 genera across all samples are given in the table below
df <- psmelt(GenusLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
df2 <- psmelt(GenusAll)
df2$Abundance<-df2$Abundance*100
GenusAll<-ddply(df2, c("Genus"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
GenusSorted<-GenusAll[order(-GenusAll$mean),]
GenusSorted[1:5,]## Genus N mean sd se
## 131 Rhodobacter 26 4.178205 4.430901 0.8689712
## 101 Methylotenera 26 2.488462 2.729264 0.5352528
## 88 Leadbetterella 26 1.992308 5.600399 1.0983287
## 58 Dysgonomonas 26 1.851282 3.369883 0.6608884
## 40 Clostridium 52 1.410897 3.143557 0.4359329GenusAll2<-ddply(df2, c("Genus","FFG"), summarise, #For Rel abu tabs
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
GenusSorted2<-GenusAll2[order(-GenusAll$mean),]
dfStats<-subset(df,FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")Trtdata2 <- ddply(df, c("Genus", "FFG","Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station)
cdataplot2
#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Genus))
#SigList<-length(unique(Trtdata$Genus))
#SigLetters2<-vector(length=NComparisons)
Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Genus", p.adjust.method = "fdr")
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Genus", p.adjust.method = "fdr")
#vec<-unlist(lst)
#for (i in levels(Means$Genus)){
# Tax<-i
# TaxAbundance<-subset(Means,Genus==i )
#print(TaxAbundance)
# Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
# difference<-TaxAbundance$p.adj
# names(difference)<-Hyphenated
# Letters<-multcompLetters(difference)
#print(Letters)
# SigList[i]<-Letters
#}
#vec<-unlist(SigList)
#vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)#+geom_text(aes(x=FFG, y=mean+se+10,label=vec))
cdataplot#+ scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(geom_text(aes(x=FFG, y=mean+se+5,label=vec))
#+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
compare_means(Abundance~FFG, data=dfStats,group.by="Genus",method = "kruskal.test",p.adjust.method="fdr")## # A tibble: 3 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacter Abundance 0.958 0.96 0.95776 ns Kruskal-Wallis
## 2 Clostridium Abundance 0.000239 0.00055 0.00024 *** Kruskal-Wallis
## 3 Methylotenera Abundance 0.000369 0.00055 0.00037 *** Kruskal-WallisMeans=compare_means(Abundance ~ FFG, data = dfStats,
group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 3 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobac~ Abunda~ Scraper Shred~ 1.00e+0 1 1.00000 ns Wilco~
## 2 Clostrid~ Abunda~ Scraper Shred~ 2.99e-4 0.00068 0.00030 *** Wilco~
## 3 Methylot~ Abunda~ Scraper Shred~ 4.57e-4 0.00068 0.00046 *** Wilco~MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats,
group.by = "Genus", p.adjust.method = "fdr")
MeansStation## # A tibble: 3 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacter Abundan~ Ostana PDR 0.0135 0.041 0.014 * Wilcox~
## 2 Clostridium Abundan~ Ostana PDR 0.325 0.49 0.325 ns Wilcox~
## 3 Methyloten~ Abundan~ Ostana PDR 0.866 0.87 0.866 ns Wilcox~# #head(Means)
# keeps <- c("Genus","group1","group2","p.adj","method","p.signif")
# keeps=Means[keeps]
# #keeps
#
#
# test3 <- list('Genus'= keeps$Genus,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
# test3= as.data.frame(test3)
# #test3
# FilteredResults<-test3[!(test3$p.adj>0.05),]
# FilteredResultsPhylumAbuAll
Below is a list of the relative abundance of all taxa at the Phylum level by FFG
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)compare_means(Abundance~FFG, data=df,group.by="Phylum",method = "kruskal.test",p.adjust.method="fdr")## # A tibble: 7 x 7
## Phylum .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Proteobacteria Abundance 0.0147 0.021 0.01466 * Kruskal-Wall~
## 2 Tenericutes Abundance 0.0643 0.075 0.06432 ns Kruskal-Wall~
## 3 Bacteroidetes Abundance 0.000400 0.0014 0.00040 *** Kruskal-Wall~
## 4 Firmicutes Abundance 0.000209 0.0014 0.00021 *** Kruskal-Wall~
## 5 Planctomycetes Abundance 0.0127 0.021 0.01274 * Kruskal-Wall~
## 6 Verrucomicrobia Abundance 0.00184 0.0043 0.00184 ** Kruskal-Wall~
## 7 Actinobacteria Abundance 0.118 0.12 0.11762 ns Kruskal-Wall~Means=compare_means(Abundance ~ FFG, data = df,
group.by = "Phylum", p.adjust.method = "fdr")Multiple comparisons for all taxa grouped by Phylum
NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Phylum))
SigList<-length(unique(Trtdata$Phylum))
SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
#for (i in levels(Means$Phylum)){
# Tax<-i
# TaxAbundance<-subset(Means,Phylum==i )
# print(TaxAbundance)
# Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
# difference<-TaxAbundance$p.adj
# names(difference)<-Hyphenated
# Letters<-multcompLetters(difference)
# print(Letters)
#
#}FamilyAbuAll
Below is a list of the relative abundance of all taxa at the Family level by FFG
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)compare_means(Abundance~FFG, data=df,group.by="Family",method = "kruskal.test",p.adjust.method="fdr")## # A tibble: 10 x 7
## Family .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobacteria~ Abundan~ 0.0518 5.20e-2 0.05176 ns Kruskal-Wa~
## 2 Aeromonadaceae Abundan~ 0.0365 4.60e-2 0.03654 * Kruskal-Wa~
## 3 Cytophagaceae Abundan~ 0.0000527 1.50e-4 5.3e-05 **** Kruskal-Wa~
## 4 Comamonadaceae Abundan~ 0.0000611 1.50e-4 6.1e-05 **** Kruskal-Wa~
## 5 Rhodobacterace~ Abundan~ 0.0464 5.20e-2 0.04640 * Kruskal-Wa~
## 6 Ruminococcaceae Abundan~ 0.0000556 1.50e-4 5.6e-05 **** Kruskal-Wa~
## 7 Desulfovibrion~ Abundan~ 0.0000686 1.50e-4 6.9e-05 **** Kruskal-Wa~
## 8 Lachnospiraceae Abundan~ 0.0000763 1.50e-4 7.6e-05 **** Kruskal-Wa~
## 9 Rikenellaceae Abundan~ 0.000129 2.10e-4 0.00013 *** Kruskal-Wa~
## 10 Methylophilace~ Abundan~ 0.000339 4.80e-4 0.00034 *** Kruskal-Wa~Means=compare_means(Abundance ~ FFG, data = df,
group.by = "Family", p.adjust.method = "fdr")Multiple comparisons for all taxa grouped by Family
NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Family))
SigList<-length(unique(Trtdata$Family))
SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
for (i in levels(Means$Family)){
Tax<-i
TaxAbundance<-subset(Means,Family==i )
print(TaxAbundance)
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
print(Letters)
}## # A tibble: 3 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Aeromonad~ Abunda~ Filter~ Scraper 0.158 0.22 0.15751 ns Wilcox~
## 2 Aeromonad~ Abunda~ Predat~ Scraper 0.0820 0.13 0.08197 ns Wilcox~
## 3 Aeromonad~ Abunda~ Scraper Shredd~ 0.0605 0.097 0.06048 ns Wilcox~
## $Letters
## Filterer Predator Scraper Shredder
## "a" "a" "a" "a"
##
## $LetterMatrix
## a
## Filterer TRUE
## Predator TRUE
## Scraper TRUE
## Shredder TRUE
##
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Comamona~ Abunda~ Filter~ Predat~ 2.40e-1 0.31 0.23953 ns Wilco~
## 2 Comamona~ Abunda~ Filter~ Scraper 6.55e-3 0.016 0.00655 ** Wilco~
## 3 Comamona~ Abunda~ Filter~ Shredd~ 1.07e-2 0.023 0.01073 * Wilco~
## 4 Comamona~ Abunda~ Predat~ Scraper 1.69e-3 0.0073 0.00169 ** Wilco~
## 5 Comamona~ Abunda~ Predat~ Shredd~ 5.28e-3 0.014 0.00528 ** Wilco~
## 6 Comamona~ Abunda~ Scraper Shredd~ 9.89e-4 0.0046 0.00099 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "a" "b" "c"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Cytophag~ Abunda~ Filter~ Predat~ 1.42e-2 0.028 0.01421 * Wilco~
## 2 Cytophag~ Abunda~ Filter~ Scraper 2.08e-3 0.0073 0.00208 ** Wilco~
## 3 Cytophag~ Abunda~ Filter~ Shredd~ 9.25e-1 0.94 0.92472 ns Wilco~
## 4 Cytophag~ Abunda~ Predat~ Scraper 7.06e-4 0.0046 0.00071 *** Wilco~
## 5 Cytophag~ Abunda~ Predat~ Shredd~ 3.41e-3 0.0095 0.00341 ** Wilco~
## 6 Cytophag~ Abunda~ Scraper Shredd~ 4.57e-4 0.0043 0.00046 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "c" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Desulfovi~ Abunda~ Filter~ Preda~ 8.91e-3 0.02 0.00891 ** Wilco~
## 2 Desulfovi~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097 *** Wilco~
## 3 Desulfovi~ Abunda~ Filter~ Shred~ 6.37e-1 0.7 0.63660 ns Wilco~
## 4 Desulfovi~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34 0.27630 ns Wilco~
## 5 Desulfovi~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259 ** Wilco~
## 6 Desulfovi~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Enterobact~ Abunda~ Filter~ Predat~ 0.0304 0.052 0.03038 * Wilco~
## 2 Enterobact~ Abunda~ Filter~ Scraper 0.0189 0.035 0.01892 * Wilco~
## 3 Enterobact~ Abunda~ Filter~ Shredd~ 0.0334 0.055 0.03336 * Wilco~
## 4 Enterobact~ Abunda~ Predat~ Scraper 0.768 0.83 0.76828 ns Wilco~
## 5 Enterobact~ Abunda~ Predat~ Shredd~ 0.886 0.94 0.88625 ns Wilco~
## 6 Enterobact~ Abunda~ Scraper Shredd~ 0.368 0.44 0.36791 ns Wilco~
## Filterer Predator Scraper Shredder
## "a" "ab" "b" "ab"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Lachnospi~ Abunda~ Filter~ Preda~ 1.21e-2 0.025 0.01206 * Wilco~
## 2 Lachnospi~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097 *** Wilco~
## 3 Lachnospi~ Abunda~ Filter~ Shred~ 6.37e-1 0.7 0.63660 ns Wilco~
## 4 Lachnospi~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34 0.27630 ns Wilco~
## 5 Lachnospi~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259 ** Wilco~
## 6 Lachnospi~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Methyloph~ Abunda~ Filter~ Preda~ 9.15e-1 0.94 0.91511 ns Wilco~
## 2 Methyloph~ Abunda~ Filter~ Scrap~ 2.08e-3 0.0073 0.00208 ** Wilco~
## 3 Methyloph~ Abunda~ Filter~ Shred~ 1.56e-1 0.22 0.15638 ns Wilco~
## 4 Methyloph~ Abunda~ Predat~ Scrap~ 8.78e-4 0.0046 0.00088 *** Wilco~
## 5 Methyloph~ Abunda~ Predat~ Shred~ 2.25e-1 0.3 0.22464 ns Wilco~
## 6 Methyloph~ Abunda~ Scraper Shred~ 4.57e-4 0.0043 0.00046 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "a" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rhodobacte~ Abunda~ Filter~ Predat~ 0.110 0.17 0.10982 ns Wilco~
## 2 Rhodobacte~ Abunda~ Filter~ Scraper 0.589 0.67 0.58915 ns Wilco~
## 3 Rhodobacte~ Abunda~ Filter~ Shredd~ 0.219 0.3 0.21930 ns Wilco~
## 4 Rhodobacte~ Abunda~ Predat~ Scraper 0.0292 0.051 0.02924 * Wilco~
## 5 Rhodobacte~ Abunda~ Predat~ Shredd~ 0.0184 0.035 0.01842 * Wilco~
## 6 Rhodobacte~ Abunda~ Scraper Shredd~ 1 1 1.00000 ns Wilco~
## Filterer Predator Scraper Shredder
## "ab" "a" "ab" "b"
## # A tibble: 5 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Rikenell~ Abunda~ Filter~ Predat~ 2.58e-2 0.047 0.02577 * Wilco~
## 2 Rikenell~ Abunda~ Filter~ Scraper 6.67e-3 0.016 0.00667 ** Wilco~
## 3 Rikenell~ Abunda~ Filter~ Shredd~ 3.95e-1 0.46 0.39509 ns Wilco~
## 4 Rikenell~ Abunda~ Predat~ Shredd~ 2.07e-3 0.0073 0.00207 ** Wilco~
## 5 Rikenell~ Abunda~ Scraper Shredd~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"
## # A tibble: 6 x 9
## Family .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Ruminococ~ Abunda~ Filter~ Preda~ 8.91e-3 0.02 0.00891 ** Wilco~
## 2 Ruminococ~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097 *** Wilco~
## 3 Ruminococ~ Abunda~ Filter~ Shred~ 1.56e-1 0.22 0.15638 ns Wilco~
## 4 Ruminococ~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34 0.27630 ns Wilco~
## 5 Ruminococ~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259 ** Wilco~
## 6 Ruminococ~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030 *** Wilco~
## Filterer Predator Scraper Shredder
## "a" "b" "b" "a"GenusAbuAll
Below is a list of the relative abundance of all taxa at the Genus level by Decomp Stage
GenusSorted2## Genus FFG N mean sd se
## 131 Caldicoprobacter Scraper 9 0.000000000 0.000000000 0.000000000
## 101 Bacteroides Filterer 4 0.000000000 0.000000000 0.000000000
## 88 Asticcacaulis Shredder 7 0.000000000 0.000000000 0.000000000
## 58 Anaerospora Predator 6 0.133333333 0.183787317 0.075030858
## 40 Agrobacterium Shredder 7 0.004761905 0.012598816 0.004761905
## 129 Caldicoprobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 133 Caulobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 69 Aquabacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 117 Bradyrhizobium Filterer 4 0.000000000 0.000000000 0.000000000
## 100 Bacillus Shredder 7 0.009523810 0.025197632 0.009523810
## 59 Anaerospora Scraper 9 0.062962963 0.176733041 0.058911014
## 121 Brevibacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 103 Bacteroides Scraper 9 0.000000000 0.000000000 0.000000000
## 6 02d06 Predator 6 0.000000000 0.000000000 0.000000000
## 70 Aquabacterium Predator 6 0.027777778 0.068041382 0.027777778
## 53 Anaerofustis Filterer 4 0.000000000 0.000000000 0.000000000
## 96 Azospirillum Shredder 7 0.000000000 0.000000000 0.000000000
## 79 Armatimonas Scraper 9 0.000000000 0.000000000 0.000000000
## 142 Chlorobaculum Predator 6 0.011111111 0.027216553 0.011111111
## 13 Achromobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 119 Bradyrhizobium Scraper 9 0.122222222 0.259272486 0.086424162
## 77 Armatimonas Filterer 4 0.008333333 0.016666667 0.008333333
## 128 Burkholderia Shredder 7 0.000000000 0.000000000 0.000000000
## 143 Chlorobaculum Scraper 9 0.000000000 0.000000000 0.000000000
## 111 Bifidobacterium Scraper 9 0.003703704 0.011111111 0.003703704
## 125 Burkholderia Filterer 4 0.000000000 0.000000000 0.000000000
## 164 CM44 Shredder 7 0.071428571 0.073102089 0.027629992
## 145 Christensenella Filterer 4 0.033333333 0.066666667 0.033333333
## 158 Clostridium Predator 12 0.000000000 0.000000000 0.000000000
## 43 Alistipes Scraper 9 0.000000000 0.000000000 0.000000000
## 73 Arenimonas Filterer 4 0.008333333 0.016666667 0.008333333
## 27 Actinomyces Scraper 9 0.000000000 0.000000000 0.000000000
## 82 Arthrobacter Predator 6 0.033333333 0.081649658 0.033333333
## 47 Anaerococcus Scraper 9 0.003703704 0.011111111 0.003703704
## 141 Chlorobaculum Filterer 4 0.000000000 0.000000000 0.000000000
## 15 Achromobacter Scraper 9 0.029629630 0.077180244 0.025726748
## 80 Armatimonas Shredder 7 0.047619048 0.066268653 0.025047197
## 17 Acidovorax Filterer 4 0.000000000 0.000000000 0.000000000
## 11 A17 Scraper 9 0.237037037 0.711111111 0.237037037
## 3 [Clostridium] Scraper 9 0.040740741 0.122222222 0.040740741
## 98 Bacillus Predator 6 0.000000000 0.000000000 0.000000000
## 157 Clostridium Filterer 8 2.845833333 4.455634475 1.575304676
## 159 Clostridium Scraper 18 0.001851852 0.007856742 0.001851852
## 94 Azospirillum Predator 6 0.138888889 0.340206909 0.138888889
## 140 Cellvibrio Shredder 7 0.042857143 0.113389342 0.042857143
## 31 Actinoplanes Scraper 9 0.000000000 0.000000000 0.000000000
## 44 Alistipes Shredder 7 0.319047619 0.264475122 0.099962200
## 81 Arthrobacter Filterer 4 0.000000000 0.000000000 0.000000000
## 20 Acidovorax Shredder 7 0.009523810 0.025197632 0.009523810
## 155 Chthoniobacter Scraper 9 0.096296296 0.233002411 0.077667470
## 54 Anaerofustis Predator 6 0.000000000 0.000000000 0.000000000
## 75 Arenimonas Scraper 9 0.000000000 0.000000000 0.000000000
## 30 Actinoplanes Predator 6 0.005555556 0.013608276 0.005555556
## 60 Anaerospora Shredder 7 0.142857143 0.280683218 0.106088285
## 55 Anaerofustis Scraper 9 0.000000000 0.000000000 0.000000000
## 132 Caldicoprobacter Shredder 7 0.028571429 0.040499526 0.015307382
## 149 Chryseobacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 118 Bradyrhizobium Predator 6 0.005555556 0.013608276 0.005555556
## 99 Bacillus Scraper 9 0.000000000 0.000000000 0.000000000
## 108 Bdellovibrio Shredder 7 0.028571429 0.052453053 0.019825390
## 39 Agrobacterium Scraper 9 0.000000000 0.000000000 0.000000000
## 35 Afifella Scraper 9 0.000000000 0.000000000 0.000000000
## 42 Alistipes Predator 6 0.000000000 0.000000000 0.000000000
## 112 Bifidobacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 85 Asticcacaulis Filterer 4 0.000000000 0.000000000 0.000000000
## 24 Acinetobacter Shredder 7 0.047619048 0.063412649 0.023967728
## 29 Actinoplanes Filterer 4 0.008333333 0.016666667 0.008333333
## 151 Chryseobacterium Scraper 9 0.000000000 0.000000000 0.000000000
## 153 Chthoniobacter Filterer 4 0.025000000 0.031914237 0.015957118
## 137 Cellvibrio Filterer 4 0.000000000 0.000000000 0.000000000
## 90 Azospira Predator 6 0.000000000 0.000000000 0.000000000
## 144 Chlorobaculum Shredder 7 0.000000000 0.000000000 0.000000000
## 56 Anaerofustis Shredder 7 0.014285714 0.037796447 0.014285714
## 126 Burkholderia Predator 6 0.000000000 0.000000000 0.000000000
## 91 Azospira Scraper 9 0.007407407 0.022222222 0.007407407
## 49 Anaerofilum Filterer 4 1.025000000 0.412198262 0.206099131
## 116 Brachybacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 154 Chthoniobacter Predator 6 0.000000000 0.000000000 0.000000000
## 41 Alistipes Filterer 4 0.000000000 0.000000000 0.000000000
## 61 Anaerotruncus Filterer 4 0.000000000 0.000000000 0.000000000
## 114 Brachybacterium Predator 6 0.000000000 0.000000000 0.000000000
## 115 Brachybacterium Scraper 9 0.092592593 0.277777778 0.092592593
## 66 Anaerovorax Predator 6 0.000000000 0.000000000 0.000000000
## 124 Brevibacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 113 Brachybacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 163 CM44 Scraper 9 0.000000000 0.000000000 0.000000000
## 19 Acidovorax Scraper 9 0.000000000 0.000000000 0.000000000
## 92 Azospira Shredder 7 0.000000000 0.000000000 0.000000000
## 161 CM44 Filterer 4 0.000000000 0.000000000 0.000000000
## 21 Acinetobacter Filterer 4 0.383333333 0.589726867 0.294863434
## 38 Agrobacterium Predator 6 0.000000000 0.000000000 0.000000000
## 62 Anaerotruncus Predator 6 0.000000000 0.000000000 0.000000000
## 120 Bradyrhizobium Shredder 7 0.000000000 0.000000000 0.000000000
## 71 Aquabacterium Scraper 9 0.000000000 0.000000000 0.000000000
## 105 Bdellovibrio Filterer 4 0.158333333 0.142400062 0.071200031
## 150 Chryseobacterium Predator 6 0.027777778 0.068041382 0.027777778
## 1 [Clostridium] Filterer 4 0.000000000 0.000000000 0.000000000
## 86 Asticcacaulis Predator 6 0.011111111 0.027216553 0.011111111
## 34 Afifella Predator 6 0.016666667 0.040824829 0.016666667
## 46 Anaerococcus Predator 6 0.000000000 0.000000000 0.000000000
## 109 Bifidobacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 162 CM44 Predator 6 0.005555556 0.013608276 0.005555556
## 16 Achromobacter Shredder 7 0.000000000 0.000000000 0.000000000
## 68 Anaerovorax Shredder 7 0.342857143 0.359379241 0.135832586
## 89 Azospira Filterer 4 0.000000000 0.000000000 0.000000000
## 102 Bacteroides Predator 6 0.000000000 0.000000000 0.000000000
## 123 Brevibacterium Scraper 9 0.129629630 0.388888889 0.129629630
## 4 [Clostridium] Shredder 7 0.000000000 0.000000000 0.000000000
## 32 Actinoplanes Shredder 7 0.000000000 0.000000000 0.000000000
## 48 Anaerococcus Shredder 7 0.000000000 0.000000000 0.000000000
## 33 Afifella Filterer 4 0.000000000 0.000000000 0.000000000
## 134 Caulobacter Predator 6 0.000000000 0.000000000 0.000000000
## 104 Bacteroides Shredder 7 0.004761905 0.012598816 0.004761905
## 18 Acidovorax Predator 6 0.000000000 0.000000000 0.000000000
## 52 Anaerofilum Shredder 7 0.519047619 0.304811504 0.115207919
## 130 Caldicoprobacter Predator 6 0.000000000 0.000000000 0.000000000
## 95 Azospirillum Scraper 9 0.000000000 0.000000000 0.000000000
## 37 Agrobacterium Filterer 4 0.000000000 0.000000000 0.000000000
## 74 Arenimonas Predator 6 0.050000000 0.078173596 0.031914237
## 146 Christensenella Predator 6 0.000000000 0.000000000 0.000000000
## 156 Chthoniobacter Shredder 7 0.004761905 0.012598816 0.004761905
## 9 A17 Filterer 4 0.000000000 0.000000000 0.000000000
## 14 Achromobacter Predator 6 0.000000000 0.000000000 0.000000000
## 67 Anaerovorax Scraper 9 0.000000000 0.000000000 0.000000000
## 87 Asticcacaulis Scraper 9 0.000000000 0.000000000 0.000000000
## 136 Caulobacter Shredder 7 0.047619048 0.125988158 0.047619048
## 160 Clostridium Shredder 14 3.611904762 4.123947307 1.102171279
## 51 Anaerofilum Scraper 9 0.000000000 0.000000000 0.000000000
## 110 Bifidobacterium Predator 6 0.000000000 0.000000000 0.000000000
## 139 Cellvibrio Scraper 9 0.000000000 0.000000000 0.000000000
## 2 [Clostridium] Predator 6 0.000000000 0.000000000 0.000000000
## 5 02d06 Filterer 4 0.000000000 0.000000000 0.000000000
## 7 02d06 Scraper 9 0.000000000 0.000000000 0.000000000
## 8 02d06 Shredder 7 0.009523810 0.025197632 0.009523810
## 22 Acinetobacter Predator 6 1.933333333 1.262097020 0.515248951
## 23 Acinetobacter Scraper 9 0.177777778 0.248886409 0.082962136
## 25 Actinomyces Filterer 4 0.000000000 0.000000000 0.000000000
## 36 Afifella Shredder 7 0.000000000 0.000000000 0.000000000
## 83 Arthrobacter Scraper 9 0.022222222 0.066666667 0.022222222
## 84 Arthrobacter Shredder 7 0.000000000 0.000000000 0.000000000
## 93 Azospirillum Filterer 4 0.000000000 0.000000000 0.000000000
## 106 Bdellovibrio Predator 6 0.244444444 0.223772371 0.091354688
## 107 Bdellovibrio Scraper 9 0.114814815 0.233399462 0.077799821
## 127 Burkholderia Scraper 9 0.029629630 0.088888889 0.029629630
## 147 Christensenella Scraper 9 0.000000000 0.000000000 0.000000000
## 10 A17 Predator 6 0.000000000 0.000000000 0.000000000
## 12 A17 Shredder 7 0.000000000 0.000000000 0.000000000
## 26 Actinomyces Predator 6 0.000000000 0.000000000 0.000000000
## 28 Actinomyces Shredder 7 0.009523810 0.025197632 0.009523810
## 45 Anaerococcus Filterer 4 0.000000000 0.000000000 0.000000000
## 50 Anaerofilum Predator 6 0.000000000 0.000000000 0.000000000
## 57 Anaerospora Filterer 4 0.033333333 0.047140452 0.023570226
## 63 Anaerotruncus Scraper 9 0.000000000 0.000000000 0.000000000
## 64 Anaerotruncus Shredder 7 0.042857143 0.065868235 0.024895853
## 65 Anaerovorax Filterer 4 0.000000000 0.000000000 0.000000000
## 72 Aquabacterium Shredder 7 0.000000000 0.000000000 0.000000000
## 76 Arenimonas Shredder 7 0.014285714 0.037796447 0.014285714
## 78 Armatimonas Predator 6 0.144444444 0.162845075 0.066481224
## 97 Bacillus Filterer 4 0.000000000 0.000000000 0.000000000
## 122 Brevibacterium Predator 6 0.038888889 0.095257934 0.038888889
## 135 Caulobacter Scraper 9 0.000000000 0.000000000 0.000000000
## 138 Cellvibrio Predator 6 0.111111111 0.226732017 0.092562958
## 148 Christensenella Shredder 7 0.000000000 0.000000000 0.000000000
## 152 Chryseobacterium Shredder 7 0.033333333 0.074535599 0.028171808PCoA
Weighted Unifrac
-Ellipses represent 95% CI for the mean of each group
physeqStats<-subset_samples(physeq, FFG!="Predator")
physeqStats<-subset_samples(physeqStats,FFG!="Filterer")
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
Jaccard
-Ellipses represent 95% CI for the mean of each group
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station,scales = "free")#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
Bray-curtis
-Ellipses represent 95% CI for the mean of each group
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
PERMANOVAs
Weighted Unifrac
Homogenieity of Multivariate Dispersions
physeqStats<-subset_samples(physeq, FFG!="Predator")
physeqStats<-subset_samples(physeqStats,FFG!="Filterer")
GPdist=phyloseq::distance(physeqStats, "wunifrac")
beta=betadisper(GPdist, sample_data(physeqStats)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00067157 0.00067157 14.985 999 0.001 ***
## Residuals 14 0.00062744 0.00004482
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
FFG
GPdist=phyloseq::distance(physeqStats, "wunifrac")
#MONMDS= ordinate(physeq, "NMDS",GPdist)adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.010177 0.0101767 4.9044 0.23867 0.002 **
## Sampling_station 1 0.003762 0.0037616 1.8128 0.08822 0.082 .
## FFG:Sampling_station 1 0.003800 0.0037998 1.8312 0.08912 0.074 .
## Residuals 12 0.024900 0.0020750 0.58399
## Total 15 0.042638 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.010177 0.0101767 4.389 0.23867 0.001 ***
## Residuals 14 0.032462 0.0023187 0.76133
## Total 15 0.042638 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Jaccard
Homogenieity of Multivariate Dispersions
GPdist=phyloseq::distance(physeq, "jaccard")
beta=betadisper(GPdist, sample_data(physeq)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.048886 0.016295 3.9989 999 0.019 *
## Residuals 22 0.089649 0.004075
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
print("Beta dispersion location")## [1] "Beta dispersion location"beta=betadisper(GPdist, sample_data(physeq)$Sampling_station)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0017764 0.00177642 2.0475 999 0.161
## Residuals 24 0.0208227 0.00086761boxplot(beta)
print("Beta dispersion Shredder v Scraper")## [1] "Beta dispersion Shredder v Scraper"physeqShredderVScraper<-subset_samples(physeq, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqShredderVScraper, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVScraper)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.010007 0.010007 3.9245 999 0.071 .
## Residuals 14 0.035700 0.002550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
print("Beta dispersion Shredder v Predator")## [1] "Beta dispersion Shredder v Predator"physeqShredderVPredator<-subset_samples(physeq, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqShredderVPredator, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVPredator)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0017264 0.0017264 0.9972 999 0.343
## Residuals 11 0.0190437 0.0017312boxplot(beta)
print("Beta dispersion Shredder v Filterer")## [1] "Beta dispersion Shredder v Filterer"physeqShredderVFilterer<-subset_samples(physeq, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqShredderVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVFilterer)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.016325 0.0163249 2.8041 999 0.155
## Residuals 9 0.052396 0.0058218boxplot(beta)
print("Beta dispersion Predator v Scraper")## [1] "Beta dispersion Predator v Scraper"physeqPredatorVScraper<-subset_samples(physeq, FFG=="Predator"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPredatorVScraper, "jaccard")
beta=betadisper(GPdist, sample_data(physeqPredatorVScraper)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.002683 0.0026825 0.936 999 0.347
## Residuals 13 0.037257 0.0028659boxplot(beta)
print("Beta dispersion Predator v Filterer")## [1] "Beta dispersion Predator v Filterer"physeqPredatorVFilterer<-subset_samples(physeq, FFG=="Predator"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqPredatorVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqPredatorVFilterer)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.025561 0.0255605 3.7911 999 0.113
## Residuals 8 0.053937 0.0067422boxplot(beta)
print("Beta dispersion Scraper v Filterer")## [1] "Beta dispersion Scraper v Filterer"physeqScraperVFilterer<-subset_samples(physeq, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqScraperVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqScraperVFilterer)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.047159 0.047159 7.3469 999 0.021 *
## Residuals 11 0.070608 0.006419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
FFG
GPdist=phyloseq::distance(physeqStats, "jaccard")
#MONMDS= ordinate(physeq, "NMDS",GPdist)adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.2666 1.26661 3.7600 0.19102 0.001 ***
## Sampling_station 1 0.6445 0.64455 1.9133 0.09721 0.012 *
## FFG:Sampling_station 1 0.6772 0.67721 2.0103 0.10213 0.013 *
## Residuals 12 4.0424 0.33687 0.60964
## Total 15 6.6308 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.2666 1.26661 3.3058 0.19102 0.001 ***
## Residuals 14 5.3642 0.38315 0.80898
## Total 15 6.6308 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Bray-curtis
Homogenieity of Multivariate Dispersions
GPdist=phyloseq::distance(physeqStats, "bray")
beta=betadisper(GPdist, sample_data(physeqStats)$FFG)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.025993 0.025993 3.8617 999 0.065 .
## Residuals 14 0.094234 0.006731
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(beta)
FFG
GPdist=phyloseq::distance(physeqStats, "bray")
#MONMDS= ordinate(physeq, "NMDS",GPdist)adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.6644 1.66440 6.4304 0.27136 0.001 ***
## Sampling_station 1 0.6540 0.65397 2.5266 0.10662 0.007 **
## FFG:Sampling_station 1 0.7092 0.70923 2.7401 0.11563 0.008 **
## Residuals 12 3.1060 0.25884 0.50639
## Total 15 6.1336 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.6644 1.66440 5.2138 0.27136 0.001 ***
## Residuals 14 4.4692 0.31923 0.72864
## Total 15 6.1336 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Jaccard Split by Sampling Station
physeqOstana<-subset_samples(physeq, Sampling_station=="Ostana")
physeqOstanaStats<-subset_samples(physeqOstana,FFG!="Predator")
physeqPDRStats<-subset_samples(physeq, Sampling_station!="Ostana")
#MONMDS= ordinate(physeq, "NMDS",GPdist)Ostana (Predator not included in stats due to N = 2)
GPdist=phyloseq::distance(physeqOstanaStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 1.8065 0.90327 2.6038 0.36653 0.001 ***
## Residuals 9 3.1222 0.34691 0.63347
## Total 11 4.9287 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper Ostana")## [1] "Shredder vs Scraper Ostana"physeqOstanaShredderVScraper<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqOstanaShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.86139 0.86139 2.3315 0.27984 0.027 *
## Residuals 6 2.21670 0.36945 0.72016
## Total 7 3.07809 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Filterer Ostana")## [1] "Shredder vs Filterer Ostana"physeqOstanaShredderVFilterer<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))## Set of permutations < 'minperm'. Generating entire set.##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.85695 0.85695 2.8363 0.36194 0.027 *
## Residuals 5 1.51071 0.30214 0.63806
## Total 6 2.36766 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Filterer Ostana")## [1] "Scraper vs Filterer Ostana"physeqOstanaScraperVFilterer<-subset_samples(physeqOstanaStats, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.9765 0.97652 2.7159 0.27953 0.011 *
## Residuals 7 2.5169 0.35956 0.72047
## Total 8 3.4934 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PDR (Includes Predator, Shredder, Scraper)
GPdist=phyloseq::distance(physeqPDRStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 1.8490 0.92448 2.8964 0.3916 0.001 ***
## Residuals 9 2.8726 0.31918 0.6084
## Total 11 4.7216 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper PDR")## [1] "Shredder vs Scraper PDR"physeqPDRShredderVScraper<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPDRShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 1.0689 1.06893 3.5129 0.36928 0.029 *
## Residuals 6 1.8257 0.30428 0.63072
## Total 7 2.8946 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Predator PDR")## [1] "Shredder vs Predator PDR"physeqPDRShredderVPredator<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.90936 0.90936 2.874 0.32387 0.022 *
## Residuals 6 1.89844 0.31641 0.67613
## Total 7 2.80779 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Predator PDR")## [1] "Scraper vs Predator PDR"physeqPDRScraperVPredator<-subset_samples(physeqPDRStats, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRScraperVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRScraperVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.79514 0.79514 2.3605 0.28234 0.02 *
## Residuals 6 2.02107 0.33684 0.71766
## Total 7 2.81621 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Random Forest
By Sampling Station
set.seed(200)
GPr = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
GPr = filter_taxa(GPr, function(x) mean(x) > 1e-5, TRUE)
PhylumAll=tax_glom(GPr, "Phylum")
FamilyAll=tax_glom(GPr,"Family")
GenusAll=tax_glom(GPr,"Genus")
ForestData=GenusAll#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$Sampling_station)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 3.85%
## Confusion matrix:
## Ostana PDR class.error
## Ostana 14 0 0.00000000
## PDR 1 11 0.08333333imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:23, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by Sampling Station")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] de46410fbc643aeaed03d6aa3878d71a X1e994910b44683e96c37da4cd04862a4
## [3] e1a9430f9c3611fc149a9f8e65bf5259 X8766563e48c605e235f4ea2485b02ee3
## [5] X768d500dec54c80d5ef769902dadd792 e6cad09c4fd44cdb1b919a77f3d92429
## [7] X237f7729fcc16a8a864b826f36f307c3 X7b9668ebf6830386ee80c560f7b275a7
## [9] a046edb519c4e3cbdf77ada497c4d743 X5e419280ad42e77b964625f286da0f08
## [11] e32f997a31db6624fc4f67c269121d4c X9eab7a1249300ee6e9757c224201f94d
## [13] b4fda0049757d0552cbf3080550ce23b X2a64f7117c72e249a9c2e5850fd6eb3e
## [15] X147463ced695e1a1dc7c68b5c64c89ed X592e7630b2e58f0bcb2a3f9cccb07ceb
## [17] X887b10e0bdad8789f41fb3c687174051 b672de3b577f55a4350eabdce0f18904
## [19] a1be228b7d7fe14508a59d112f71fa1a X129c4eff55474784a99b0c7309c15c50
## [21] ca23f159cd2907c86ea8d9cf10ae5df9 X0fbb434b42e7241efa6451d8f1b429d0
## [23] ad07e5885874179952e16bf2f7bb57b3
## 165 Levels: de46410fbc643aeaed03d6aa3878d71a ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| b4fda0049757d0552cbf3080550ce23b | Bacteria | Proteobacteria | Alphaproteobacteria | Rhodobacterales | Rhodobacteraceae | Anaerospora | NA | NA | NA | NA | NA | NA | NA | NA |
| e1a9430f9c3611fc149a9f8e65bf5259 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhodobacterales | Rhodobacteraceae | Rhodobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| ca23f159cd2907c86ea8d9cf10ae5df9 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Hyphomicrobiaceae | Hyphomicrobium | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| e6cad09c4fd44cdb1b919a77f3d92429 | Bacteria | Bacteroidetes | [Saprospirae] | [Saprospirales] | Chitinophagaceae | Sediminibacterium | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
| a1be228b7d7fe14508a59d112f71fa1a | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Rudanella | NA | NA | NA | NA | NA | NA | NA | NA |
| de46410fbc643aeaed03d6aa3878d71a | Bacteria | Verrucomicrobia | Verrucomicrobiae | Verrucomicrobiales | Verrucomicrobiaceae | Luteolibacter | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "Sampling_station"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
Trtdata## Genus Sampling_station N mean sd se
## 1 Anaerospora Ostana 14 0.16666667 0.24564584 0.065651613
## 2 Anaerospora PDR 12 0.01388889 0.03320683 0.009585986
## 3 Emticicia Ostana 14 0.50952381 0.54685247 0.146152469
## 4 Emticicia PDR 12 1.06666667 1.39544714 0.402830892
## 5 Hyphomicrobium Ostana 14 0.33571429 0.45750425 0.122273153
## 6 Hyphomicrobium PDR 12 0.18333333 0.40489430 0.116882916
## 7 Leadbetterella Ostana 14 0.14523810 0.26718128 0.071407201
## 8 Leadbetterella PDR 12 4.14722222 7.86069956 2.269188505
## 9 Luteolibacter Ostana 14 0.41904762 0.48739922 0.130262920
## 10 Luteolibacter PDR 12 0.28333333 0.45979794 0.132732231
## 11 Perlucidibaca Ostana 14 0.34047619 0.45726399 0.122208941
## 12 Perlucidibaca PDR 12 1.93055556 3.20673087 0.925703467
## 13 Rhodobacter Ostana 14 4.97142857 5.19945522 1.389612859
## 14 Rhodobacter PDR 12 3.25277778 3.30687861 0.954613627
## 15 Rhodoferax Ostana 14 1.14761905 1.12376917 0.300339943
## 16 Rhodoferax PDR 12 1.31666667 1.24182173 0.358483055
## 17 Rudanella Ostana 14 0.16428571 0.36408515 0.097305850
## 18 Rudanella PDR 12 0.01944444 0.06735753 0.019444444
## 19 Sediminibacterium Ostana 14 0.05476190 0.08534209 0.022808633
## 20 Sediminibacterium PDR 12 0.26944444 0.45736583 0.132030142cdataplot=ggplot(Trtdata, aes(x=Sampling_station,y=mean))+geom_bar(aes(fill = Sampling_station),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Sampling Station")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
compare_means(Abundance ~ Sampling_station, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterella Abundance 0.250 0.5 0.250 ns Kruskal-Wall~
## 2 Rhodobacter Abundance 0.136 0.34 0.136 ns Kruskal-Wall~
## 3 Perlucidibaca Abundance 0.660 0.73 0.660 ns Kruskal-Wall~
## 4 Rhodoferax Abundance 0.796 0.8 0.796 ns Kruskal-Wall~
## 5 Emticicia Abundance 0.531 0.66 0.531 ns Kruskal-Wall~
## 6 Luteolibacter Abundance 0.0435 0.15 0.044 * Kruskal-Wall~
## 7 Sediminibacterium Abundance 0.304 0.51 0.304 ns Kruskal-Wall~
## 8 Hyphomicrobium Abundance 0.411 0.59 0.411 ns Kruskal-Wall~
## 9 Rudanella Abundance 0.0350 0.15 0.035 * Kruskal-Wall~
## 10 Anaerospora Abundance 0.0276 0.15 0.028 * Kruskal-Wall~ByFFGOstana
GenusAllOstanaStats<-subset_samples(GenusAll, Sampling_station=="Ostana")
GenusAllOstanaStats<-subset_samples(GenusAllOstanaStats, FFG!="Predator")
ForestData=GenusAllOstanaStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 8.33%
## Confusion matrix:
## Filterer Scraper Shredder class.error
## Filterer 3 1 0 0.25
## Scraper 0 5 0 0.00
## Shredder 0 0 3 0.00imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by FFG")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] acccb7cec4d146864bc11d37da55dcd0 a7c725b951d62e6f6814fb8ca64a356e
## [3] ad07e5885874179952e16bf2f7bb57b3 X7403bfa23d688a94b469051d9f13ae5f
## [5] X0fbb434b42e7241efa6451d8f1b429d0 b672de3b577f55a4350eabdce0f18904
## [7] a42594aa39e8c4b1efb497c275e791ba a046edb519c4e3cbdf77ada497c4d743
## [9] X7b9668ebf6830386ee80c560f7b275a7 dc7d0d97c9604c01c99c972b878e09a0
## [11] X8766563e48c605e235f4ea2485b02ee3 X07705f7fffe33ad226eab098040fbb75
## [13] X1036caae06a5e5c605ab8120e5b5b7bc e32f997a31db6624fc4f67c269121d4c
## [15] X43dc593942c5c838cdcb10dac0a39474 d9a11b1b600813533e2f5f0cd7a2eb44
## [17] X4367180b0dc40f77a063355ebab6a4bb ca23f159cd2907c86ea8d9cf10ae5df9
## [19] X855eedb949870fe66f915b9293e0de50
## 165 Levels: acccb7cec4d146864bc11d37da55dcd0 ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| a7c725b951d62e6f6814fb8ca64a356e | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Clostridium | NA | NA | NA | NA | NA | NA | NA | NA |
| d9a11b1b600813533e2f5f0cd7a2eb44 | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Coprococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| dc7d0d97c9604c01c99c972b878e09a0 | Bacteria | Firmicutes | Clostridia | Clostridiales | Ruminococcaceae | Ruminococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| ca23f159cd2907c86ea8d9cf10ae5df9 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Hyphomicrobiaceae | Hyphomicrobium | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| acccb7cec4d146864bc11d37da55dcd0 | Bacteria | Proteobacteria | Betaproteobacteria | Methylophilales | Methylophilaceae | Methylotenera | NA | NA | NA | NA | NA | NA | NA | NA |
| a42594aa39e8c4b1efb497c275e791ba | Bacteria | Deferribacteres | Deferribacteres | Deferribacterales | Deferribacteraceae | Mucispirillum | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterel~ Abundance 0.000133 0.00027 0.00013 *** Kruskal-Wal~
## 2 Clostridium Abundance 0.0000465 0.00023 4.6e-05 **** Kruskal-Wal~
## 3 Perlucidibaca Abundance 0.0000695 0.00023 7.0e-05 **** Kruskal-Wal~
## 4 Methylotenera Abundance 0.000339 0.00048 0.00034 *** Kruskal-Wal~
## 5 Mucispirillum Abundance 0.000268 0.00045 0.00027 *** Kruskal-Wal~
## 6 Rhodoferax Abundance 0.000816 0.001 0.00082 *** Kruskal-Wal~
## 7 Emticicia Abundance 0.000131 0.00027 0.00013 *** Kruskal-Wal~
## 8 Ruminococcus Abundance 0.0000225 0.00023 2.3e-05 **** Kruskal-Wal~
## 9 Coprococcus Abundance 0.00663 0.0074 0.00663 ** Kruskal-Wal~
## 10 Hyphomicrobi~ Abundance 0.140 0.14 0.13988 ns Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 53 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.025 0.01392 * Wilco~
## 2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.015 0.00653 ** Wilco~
## 3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.51 0.44341 ns Wilco~
## 4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0047 0.00043 *** Wilco~
## 5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.013 0.00528 ** Wilco~
## 6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0052 0.00125 ** Wilco~
## 7 Clostrid~ Abunda~ Filter~ Preda~ 5.74e-3 0.014 0.00574 ** Wilco~
## 8 Clostrid~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0047 0.00097 *** Wilco~
## 9 Clostrid~ Abunda~ Filter~ Shred~ 7.77e-1 0.79 0.77681 ns Wilco~
## 10 Clostrid~ Abunda~ Predat~ Shred~ 2.07e-3 0.0061 0.00207 ** Wilco~
## # ... with 43 more rowsMeans=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
Tax<-i
TaxAbundance<-subset(Means,Genus==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%) Forest") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot
ByFFGPDR
GenusAllPDRStats<-subset_samples(GenusAll, Sampling_station!="Ostana")
ForestData=GenusAllPDRStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 0%
## Confusion matrix:
## Predator Scraper Shredder class.error
## Predator 4 0 0 0
## Scraper 0 4 0 0
## Shredder 0 0 4 0imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by FFG")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] a046edb519c4e3cbdf77ada497c4d743 X2b05a6e9e6c580f7fc168a5bbb29de2d
## [3] acccb7cec4d146864bc11d37da55dcd0 ad07e5885874179952e16bf2f7bb57b3
## [5] X237f7729fcc16a8a864b826f36f307c3 X5e419280ad42e77b964625f286da0f08
## [7] X1036caae06a5e5c605ab8120e5b5b7bc X1e994910b44683e96c37da4cd04862a4
## [9] X07705f7fffe33ad226eab098040fbb75 X6b78b6e57cad2b7f9420cee2c46db2aa
## [11] edb114398ced5bf958f4f404493a6642 e32f997a31db6624fc4f67c269121d4c
## [13] a7c725b951d62e6f6814fb8ca64a356e X4367180b0dc40f77a063355ebab6a4bb
## [15] e1a9430f9c3611fc149a9f8e65bf5259 de46410fbc643aeaed03d6aa3878d71a
## [17] b672de3b577f55a4350eabdce0f18904 dc7d0d97c9604c01c99c972b878e09a0
## [19] X0fbb434b42e7241efa6451d8f1b429d0
## 165 Levels: a046edb519c4e3cbdf77ada497c4d743 ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| a7c725b951d62e6f6814fb8ca64a356e | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Clostridium | NA | NA | NA | NA | NA | NA | NA | NA |
| dc7d0d97c9604c01c99c972b878e09a0 | Bacteria | Firmicutes | Clostridia | Clostridiales | Ruminococcaceae | Ruminococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| e1a9430f9c3611fc149a9f8e65bf5259 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhodobacterales | Rhodobacteraceae | Rhodobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| acccb7cec4d146864bc11d37da55dcd0 | Bacteria | Proteobacteria | Betaproteobacteria | Methylophilales | Methylophilaceae | Methylotenera | NA | NA | NA | NA | NA | NA | NA | NA |
| edb114398ced5bf958f4f404493a6642 | Bacteria | Bacteroidetes | [Saprospirae] | [Saprospirales] | Saprospiraceae | Haliscomenobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
| de46410fbc643aeaed03d6aa3878d71a | Bacteria | Verrucomicrobia | Verrucomicrobiae | Verrucomicrobiales | Verrucomicrobiaceae | Luteolibacter | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
#cdataplot
compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterella Abundan~ 1.33e-4 0.00022 0.00013 *** Kruskal-Wal~
## 2 Rhodobacter Abundan~ 4.63e-2 0.046 0.04634 * Kruskal-Wal~
## 3 Clostridium Abundan~ 4.65e-5 0.00022 4.6e-05 **** Kruskal-Wal~
## 4 Perlucidibaca Abundan~ 6.95e-5 0.00022 7.0e-05 **** Kruskal-Wal~
## 5 Methylotenera Abundan~ 3.39e-4 0.00048 0.00034 *** Kruskal-Wal~
## 6 Rhodoferax Abundan~ 8.16e-4 0.001 0.00082 *** Kruskal-Wal~
## 7 Emticicia Abundan~ 1.31e-4 0.00022 0.00013 *** Kruskal-Wal~
## 8 Ruminococcus Abundan~ 2.25e-5 0.00022 2.3e-05 **** Kruskal-Wal~
## 9 Luteolibacter Abundan~ 1.09e-2 0.012 0.01091 * Kruskal-Wal~
## 10 Haliscomenobac~ Abundan~ 1.16e-4 0.00022 0.00012 *** Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 55 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.025 0.01392 * Wilco~
## 2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.015 0.00653 ** Wilco~
## 3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.51 0.44341 ns Wilco~
## 4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0042 0.00043 *** Wilco~
## 5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.014 0.00528 ** Wilco~
## 6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0057 0.00125 ** Wilco~
## 7 Rhodobac~ Abunda~ Filter~ Preda~ 1.10e-1 0.15 0.10982 ns Wilco~
## 8 Rhodobac~ Abunda~ Filter~ Scrap~ 5.89e-1 0.65 0.58915 ns Wilco~
## 9 Rhodobac~ Abunda~ Filter~ Shred~ 2.18e-1 0.28 0.21825 ns Wilco~
## 10 Rhodobac~ Abunda~ Predat~ Scrap~ 2.92e-2 0.046 0.02924 * Wilco~
## # ... with 45 more rowsMeans=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
Tax<-i
TaxAbundance<-subset(Means,Genus==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%) Alpine Prairie") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot
ByFFGCombined
#GenusAllStats<-subset_samples(GenusAll, FFG!="Predator")
#GenusAllStats<-subset_samples(GenusAllStats, FFG!="Filterer")
GenusAllStats<-GenusAll
ForestData=GenusAllStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results##
## Call:
## randomForest(formula = response ~ ., data = rf.data, ntree = 1000)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 12
##
## OOB estimate of error rate: 3.85%
## Confusion matrix:
## Filterer Predator Scraper Shredder class.error
## Filterer 3 1 0 0 0.25
## Predator 0 6 0 0 0.00
## Scraper 0 0 9 0 0.00
## Shredder 0 0 0 7 0.00imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
geom_bar(stat = "identity", fill = "indianred") +
coord_flip() +
ggtitle("Most important genera for classifying samples\n by FFG")#\n in a string tells it to start a new line
#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames## [1] ad07e5885874179952e16bf2f7bb57b3 X0fbb434b42e7241efa6451d8f1b429d0
## [3] b672de3b577f55a4350eabdce0f18904 X43dc593942c5c838cdcb10dac0a39474
## [5] e32f997a31db6624fc4f67c269121d4c acccb7cec4d146864bc11d37da55dcd0
## [7] X4367180b0dc40f77a063355ebab6a4bb a7c725b951d62e6f6814fb8ca64a356e
## [9] a046edb519c4e3cbdf77ada497c4d743 dc7d0d97c9604c01c99c972b878e09a0
## [11] X2b05a6e9e6c580f7fc168a5bbb29de2d a42594aa39e8c4b1efb497c275e791ba
## [13] edb114398ced5bf958f4f404493a6642 X1036caae06a5e5c605ab8120e5b5b7bc
## [15] X7403bfa23d688a94b469051d9f13ae5f X6b78b6e57cad2b7f9420cee2c46db2aa
## [17] ef53cef3b90f7b7b0189b52eaf440bfd X237f7729fcc16a8a864b826f36f307c3
## [19] X07705f7fffe33ad226eab098040fbb75
## 165 Levels: ad07e5885874179952e16bf2f7bb57b3 ...PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable| Kingdom | Phylum | Class | Order | Family | Genus | Species | Rank5 | Rank6 | Rank7 | Rank4 | Rank2 | Rank3 | Rank1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| a7c725b951d62e6f6814fb8ca64a356e | Bacteria | Firmicutes | Clostridia | Clostridiales | Lachnospiraceae | Clostridium | NA | NA | NA | NA | NA | NA | NA | NA |
| dc7d0d97c9604c01c99c972b878e09a0 | Bacteria | Firmicutes | Clostridia | Clostridiales | Ruminococcaceae | Ruminococcus | NA | NA | NA | NA | NA | NA | NA | NA |
| b672de3b577f55a4350eabdce0f18904 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Perlucidibaca | NA | NA | NA | NA | NA | NA | NA | NA |
| a046edb519c4e3cbdf77ada497c4d743 | Bacteria | Proteobacteria | Betaproteobacteria | Burkholderiales | Comamonadaceae | Rhodoferax | NA | NA | NA | NA | NA | NA | NA | NA |
| acccb7cec4d146864bc11d37da55dcd0 | Bacteria | Proteobacteria | Betaproteobacteria | Methylophilales | Methylophilaceae | Methylotenera | NA | NA | NA | NA | NA | NA | NA | NA |
| a42594aa39e8c4b1efb497c275e791ba | Bacteria | Deferribacteres | Deferribacteres | Deferribacterales | Deferribacteraceae | Mucispirillum | NA | NA | NA | NA | NA | NA | NA | NA |
| ef53cef3b90f7b7b0189b52eaf440bfd | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Runella | NA | NA | NA | NA | NA | NA | NA | NA |
| edb114398ced5bf958f4f404493a6642 | Bacteria | Bacteroidetes | [Saprospirae] | [Saprospirales] | Saprospiraceae | Haliscomenobacter | NA | NA | NA | NA | NA | NA | NA | NA |
| e32f997a31db6624fc4f67c269121d4c | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Leadbetterella | NA | NA | NA | NA | NA | NA | NA | NA |
| ad07e5885874179952e16bf2f7bb57b3 | Bacteria | Bacteroidetes | Cytophagia | Cytophagales | Cytophagaceae | Emticicia | NA | NA | NA | NA | NA | NA | NA | NA |
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 26 samples ]
## sample_data() Sample Data: [ 26 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 14 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 10 tips and 9 internal nodes ]df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
N = length(Abundance),
mean = mean(Abundance),
sd = sd(Abundance),
se = sd / sqrt(N)
)
Trtdata## Genus FFG N mean sd se
## 1 Clostridium Filterer 4 5.68333333 4.98520032 2.49260016
## 2 Clostridium Predator 6 0.00000000 0.00000000 0.00000000
## 3 Clostridium Scraper 9 0.00000000 0.00000000 0.00000000
## 4 Clostridium Shredder 7 6.82380952 3.53619980 1.33655789
## 5 Emticicia Filterer 4 0.76666667 0.36004115 0.18002057
## 6 Emticicia Predator 6 2.06666667 1.42439071 0.58150507
## 7 Emticicia Scraper 9 0.01111111 0.03333333 0.01111111
## 8 Emticicia Shredder 7 0.62380952 0.47442530 0.17931591
## 9 Haliscomenobacter Filterer 4 0.07500000 0.07391186 0.03695593
## 10 Haliscomenobacter Predator 6 0.13888889 0.12003086 0.04900239
## 11 Haliscomenobacter Scraper 9 0.00000000 0.00000000 0.00000000
## 12 Haliscomenobacter Shredder 7 0.00000000 0.00000000 0.00000000
## 13 Leadbetterella Filterer 4 0.10000000 0.11547005 0.05773503
## 14 Leadbetterella Predator 6 8.38333333 9.64634070 3.93810210
## 15 Leadbetterella Scraper 9 0.00000000 0.00000000 0.00000000
## 16 Leadbetterella Shredder 7 0.15714286 0.14104673 0.05331065
## 17 Methylotenera Filterer 4 2.50833333 1.45178689 0.72589345
## 18 Methylotenera Predator 6 3.15555556 3.02571693 1.23524377
## 19 Methylotenera Scraper 9 0.01111111 0.03333333 0.01111111
## 20 Methylotenera Shredder 7 5.09047619 2.19145768 0.82829315
## 21 Mucispirillum Filterer 4 2.76666667 2.05065482 1.02532741
## 22 Mucispirillum Predator 6 0.01666667 0.04082483 0.01666667
## 23 Mucispirillum Scraper 9 0.00000000 0.00000000 0.00000000
## 24 Mucispirillum Shredder 7 0.81428571 0.57214902 0.21625200
## 25 Perlucidibaca Filterer 4 0.75833333 0.30107831 0.15053915
## 26 Perlucidibaca Predator 6 4.11111111 3.41107130 1.39256403
## 27 Perlucidibaca Scraper 9 0.00000000 0.00000000 0.00000000
## 28 Perlucidibaca Shredder 7 0.03333333 0.05091751 0.01924501
## 29 Rhodoferax Filterer 4 2.19166667 1.59591492 0.79795746
## 30 Rhodoferax Predator 6 2.16111111 0.97260627 0.39706485
## 31 Rhodoferax Scraper 9 0.18518519 0.32236587 0.10745529
## 32 Rhodoferax Shredder 7 1.20952381 0.52200266 0.19729846
## 33 Ruminococcus Filterer 4 0.00000000 0.00000000 0.00000000
## 34 Ruminococcus Predator 6 0.00000000 0.00000000 0.00000000
## 35 Ruminococcus Scraper 9 0.00000000 0.00000000 0.00000000
## 36 Ruminococcus Shredder 7 1.16666667 1.20308246 0.45472243
## 37 Runella Filterer 4 0.27500000 0.44253060 0.22126530
## 38 Runella Predator 6 0.70555556 0.91394546 0.37311667
## 39 Runella Scraper 9 0.00000000 0.00000000 0.00000000
## 40 Runella Shredder 7 0.01428571 0.03779645 0.01428571cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")## # A tibble: 10 x 7
## Genus .y. p p.adj p.format p.signif method
## <fct> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbetterella Abundan~ 1.33e-4 0.00022 0.00013 *** Kruskal-Wal~
## 2 Clostridium Abundan~ 4.65e-5 0.00022 4.6e-05 **** Kruskal-Wal~
## 3 Perlucidibaca Abundan~ 6.95e-5 0.00022 7.0e-05 **** Kruskal-Wal~
## 4 Methylotenera Abundan~ 3.39e-4 0.00038 0.00034 *** Kruskal-Wal~
## 5 Mucispirillum Abundan~ 2.68e-4 0.00036 0.00027 *** Kruskal-Wal~
## 6 Rhodoferax Abundan~ 8.16e-4 0.00082 0.00082 *** Kruskal-Wal~
## 7 Emticicia Abundan~ 1.31e-4 0.00022 0.00013 *** Kruskal-Wal~
## 8 Ruminococcus Abundan~ 2.25e-5 0.00022 2.3e-05 **** Kruskal-Wal~
## 9 Runella Abundan~ 2.89e-4 0.00036 0.00029 *** Kruskal-Wal~
## 10 Haliscomenobac~ Abundan~ 1.16e-4 0.00022 0.00012 *** Kruskal-Wal~Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means## # A tibble: 55 x 9
## Genus .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.022 0.01392 * Wilco~
## 2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.013 0.00653 ** Wilco~
## 3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.48 0.44341 ns Wilco~
## 4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0036 0.00043 *** Wilco~
## 5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.012 0.00528 ** Wilco~
## 6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0046 0.00125 ** Wilco~
## 7 Clostrid~ Abunda~ Filter~ Preda~ 5.74e-3 0.012 0.00574 ** Wilco~
## 8 Clostrid~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0041 0.00097 *** Wilco~
## 9 Clostrid~ Abunda~ Filter~ Shred~ 7.77e-1 0.79 0.77681 ns Wilco~
## 10 Clostrid~ Abunda~ Predat~ Shred~ 2.07e-3 0.0054 0.00207 ** Wilco~
## # ... with 45 more rowsMeans=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
Tax<-i
TaxAbundance<-subset(Means,Genus==i )
Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
difference<-TaxAbundance$p.adj
names(difference)<-Hyphenated
Letters<-multcompLetters(difference)
#print(Letters)
SigList[i]<-Letters
}## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement lengthvec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot
PiCRUST
PiCRUST<-import_biom("C:\\Users\\Joe Receveur\\Downloads\\ItalyInvertPiCRUST2.20.19.biom")
PiCRUST<-merge_phyloseq(PiCRUST,sampdat)
PiCRUST=rarefy_even_depth(PiCRUST, 750000, replace = TRUE, trimOTUs = TRUE, verbose = TRUE,rngseed = TRUE)## `set.seed(TRUE)` was used to initialize repeatable random subsampling.## Please record this for your records so others can reproduce.## Try `set.seed(TRUE); .Random.seed` for the full vector## ...## 487OTUs were removed because they are no longer
## present in any sample after random subsampling## ...PiCRUST by sampling x FFG
GPdist=phyloseq::distance(PiCRUST, "jaccard")
adonis(GPdist ~ FFG*Sampling_station, as(sample_data(PiCRUST), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(PiCRUST), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 3 0.90035 0.300117 5.0869 0.40838 0.001 ***
## Sampling_station 1 0.08310 0.083101 1.4085 0.03769 0.175
## FFG:Sampling_station 2 0.10029 0.050145 0.8499 0.04549 0.641
## Residuals 19 1.12096 0.058998 0.50844
## Total 25 2.20470 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PiCRUST beta div by FFG
print("Shredder vs Scraper")## [1] "Shredder vs Scraper"PiCRUSTShredderVScraper<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(PiCRUSTShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.50616 0.50616 6.4478 0.31533 0.001 ***
## Residuals 14 1.09902 0.07850 0.68467
## Total 15 1.60518 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Filterer")## [1] "Shredder vs Filterer"PiCRUSTShredderVFilterer<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTShredderVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTShredderVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.057191 0.057191 2.5935 0.2237 0.027 *
## Residuals 9 0.198462 0.022051 0.7763
## Total 10 0.255653 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Predator")## [1] "Shredder vs Predator"physeqShredderVPredator<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqShredderVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqShredderVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.33225 0.33225 17.331 0.61173 0.002 **
## Residuals 11 0.21088 0.01917 0.38827
## Total 12 0.54313 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Filterer")## [1] "Scraper vs Filterer"PiCRUSTScraperVFilterer<-subset_samples(PiCRUST, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTScraperVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTScraperVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.29397 0.293967 2.9572 0.21188 0.014 *
## Residuals 11 1.09347 0.099406 0.78812
## Total 12 1.38744 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Predator")## [1] "Scraper vs Predator"physeqScraperVPredator<-subset_samples(PiCRUST, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqScraperVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqScraperVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.3455 0.34550 4.0614 0.23805 0.002 **
## Residuals 13 1.1059 0.08507 0.76195
## Total 14 1.4514 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Predator vs Filterer")## [1] "Predator vs Filterer"PiCRUSTPredatorVFilterer<-subset_samples(PiCRUST, FFG=="Predator"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTPredatorVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTPredatorVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTPredatorVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.13581 0.135811 5.2914 0.39811 0.013 *
## Residuals 8 0.20533 0.025666 0.60189
## Total 9 0.34114 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PiCRUST Jaccard Split by Sampling Station
physeqOstana<-subset_samples(PiCRUST, Sampling_station=="Ostana")
physeqOstanaStats<-subset_samples(physeqOstana,FFG!="Predator")
physeqPDRStats<-subset_samples(PiCRUST, Sampling_station!="Ostana")
#MONMDS= ordinate(physeq, "NMDS",GPdist)Ostana (Predator not included in stats due to N = 2)
GPdist=phyloseq::distance(physeqOstanaStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 0.32232 0.161160 2.2914 0.3374 0.018 *
## Residuals 9 0.63299 0.070333 0.6626
## Total 11 0.95531 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper Ostana")## [1] "Shredder vs Scraper Ostana"physeqOstanaShredderVScraper<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqOstanaShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.20183 0.201826 2.257 0.27334 0.11
## Residuals 6 0.53654 0.089423 0.72666
## Total 7 0.73836 1.00000print("Shredder vs Filterer Ostana")## [1] "Shredder vs Filterer Ostana"physeqOstanaShredderVFilterer<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))## Set of permutations < 'minperm'. Generating entire set.##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.03489 0.034890 1.2951 0.20573 0.294
## Residuals 5 0.13470 0.026939 0.79427
## Total 6 0.16959 1.00000print("Scraper vs Filterer Ostana")## [1] "Scraper vs Filterer Ostana"physeqOstanaScraperVFilterer<-subset_samples(physeqOstanaStats, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.22322 0.223224 2.6273 0.2729 0.024 *
## Residuals 7 0.59475 0.084965 0.7271
## Total 8 0.81798 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PDR (Includes Predator, Shredder, Scraper)
GPdist=phyloseq::distance(physeqPDRStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRStats), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRStats), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 2 0.61094 0.305470 6.0445 0.57324 0.001 ***
## Residuals 9 0.45483 0.050537 0.42676
## Total 11 1.06577 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Scraper PDR")## [1] "Shredder vs Scraper PDR"physeqPDRShredderVScraper<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPDRShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVScraper), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVScraper), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.35923 0.35923 5.2268 0.46556 0.041 *
## Residuals 6 0.41237 0.06873 0.53444
## Total 7 0.77161 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Shredder vs Predator PDR")## [1] "Shredder vs Predator PDR"physeqPDRShredderVPredator<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.27869 0.278694 24.559 0.80366 0.04 *
## Residuals 6 0.06809 0.011348 0.19634
## Total 7 0.34678 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print("Scraper vs Predator PDR")## [1] "Scraper vs Predator PDR"physeqPDRScraperVPredator<-subset_samples(physeqPDRStats, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRScraperVPredator), "data.frame"))##
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRScraperVPredator), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## FFG 1 0.27848 0.278483 3.8931 0.39352 0.028 *
## Residuals 6 0.42920 0.071533 0.60648
## Total 7 0.70768 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PCoA PiCRUST
ord=ordinate(PiCRUST,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ord=ordinate(PiCRUST,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse